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I.         OVERVIEW  OF  THESIS 
A.      THESIS  OBJECTIVES 

The  development  of  the  Simplex  algorithm  by  Dantzig  in  middle  of  this  century 
represents  a  milestone  in  linear  optimization  techniques.  The  impact  of  Dantzig's 
work  is  profound.  Results  of  his  work  include  the  revival  or  introduction  of  a  number 
of  mathematical  disciplines,  including  convexity  and  duality  theories.  Applications 
for  the  Simplex  algorithm,  and  the  accompanying  refinements,  are  vast,  and  many 
continue  to  explore  new  and  diverse  applications. 

The  majority  of  the  research  on  linear  optimization  problems  is  taking  place 
in  various  fields  of  Operations  Research.  Of  course,  the  Simplex  algorithm  itself 
is  particularly  well  suited  to  problems  in  that  particular  discipline,  rendering  rapid 
solutions  to  production  planning  models,  transportation  problems,  and  a  variety  of 
other  "real  world"  applications.  A  great  deal  of  work  was  done  up  to  the  early 
1970's  in  attempts  to  mold  the  Simplex  algorithm  into  an  engineering  and  theoretical 
mathematical  tool.  With  the  advent  of  more  sophisticated  computer  hardware  and 
software,  there  may  be  utility  in  reconsidering  the  role  of  the  Simplex  algorithm  in 
control,  approximation,  and  other  infinite  dimensional  applications. 

This  document  is  intended  to  serve  two  main  purposes.  First,  the  thesis  is 
intended  to  serve  as  an  introduction  to  linear  optimization  and  to  the  Simplex  algo- 
rithm, or  a  theoretical  review  for  readers  already  familiar  with  these  topics.  Second, 


it  is  intended  to  present  less  traditional  problems  in  a  manner  that  is  suitable  for 
solution  with  the  Simplex  algorithm. 

B.      THESIS  FORMAT 

The  thesis  is  broken  into  three  parts.  The  first  part,  consisting  of  the  first 
two  chapters,  is  devoted  to  describing  sample  problems  with  which  the  theory  of  the 
Simplex  algorithm  is  illustrated.  Also  image  reconstruction  is  introduced,  a  problem 
whose  solution  by  the  Simplex  algorithm  highlights  the  thesis.  These  examples  are 
more  fully  developed  in  the  latter  sections. 

The  first  example  is  particularly  unusual,  as  we  find  an  orthogonal  basis  of 
the  infinite  dimensional  vector  space  L2[0, 1].  To  the  author's  knowledge,  this  is  the 
first  attempt  to  use  the  Simplex  algorithm  in  this  capacity.  The  formulations  that 
result  from  this  problem  are  particularly  easy  to  understand,  and  lend  a  great  deal 
of  understanding  to  concepts  underlying  the  Simplex  algorithm. 

The  second  example  may  be  found  infrequently  in  literature  on  linear  opti- 
mization. We  seek  the  best  approximation  to  the  exponential  function  over  a  closed 
interval  in  the  uniform  norm  sense.  That  is,  we  formulate  a  uniform  approximation 
problem  as  a  linear  optimization  problem.  The  formulation  is  used  primarily  in  the 
discussion  of  duality. 

The  final  example  is  again  a  novel  one.  We  formulate  the  problem  of  computer 
assisted  tomographic  (CAT)  image  reconstruction  as  a  linear  optimization  problem, 
and  solve  a  small  sample  problem  with  the  Simplex  algorithm. 


The  second  portion,  consisting  of  Chapters  IV  and  V,  introduces  the  machinery 
behind  the  Simplex  algorithm,  culminating  with  a  brief  introduction  to  the  algorithm 
itself.  Chapter  IV  is  an  exploration  of  convexity,  both  as  it  pertains  to  sets  and 
functions.  The  major  emphasis  of  the  chapter  is  on  convex  subsets  of  TV1 .  Chapter  V 
builds  on  the  convexity  results  as  they  pertain  to  duality.  Fundamental  concepts  of 
duality  are  presented  in  this  chapter,  and  it  concludes  with  a  generic  description  of 
the  algorithm. 

The  thesis  concludes  with  the  formulation  of  the  image  reconstruction  problem 
as  a  linear  optimization  problem  in  the  general  case.  The  first  portion  of  the  chapter 
is  devoted  to  the  formulation,  followed  by  the  statement  of  the  dual  problem.  Finally, 
a  sample  problem  is  solved,  and  some  analysis  of  the  appropriateness  of  the  Simplex 
algorithm  as  a  solution  tool  for  this  particular  problem  is  offered. 


II.         PRELIMINARIES 

A.  OVERVIEW 

We  devote  this  chapter  to  the  preliminaries  of  linear  optimization.  We  begin 
by  defining  three  very  different  examples,  which  we  develop  as  a  means  to  explore 
linear  optimization  methods.  We  then  define  the  optimization  problem  in  general, 
and  the  linear  optimization  problem  specifically.  We  close  with  a  synopsis  of  the 
assumptions  that  characterize  the  linear  optimization  problem. 

B.  FIRST  EXAMPLES 

This  thesis  extensively  discusses  three  examples.  We  begin  by  stating  two  of 
our  three  examples  to  which  we  refer  throughout  the  thesis.  Because  of  its  complexity 
and  importance  to  this  work,  the  third  example  is  treated  separately. 

1.       Example  1:  Generation  of  an  Orthogonal  Basis  for 

I2[0, 1] 

Our  first  example  is  one  of  importance  in  many  areas  of  approximation. 
We  wish  to  find  some  orthogonal  basis  for  an  infinite  dimensional  vector  space.  The 
utility  of  such  bases  may  be  found  in  any  elementary  linear  algebra  or  applied  math- 
ematics text.  The  interested  reader  is  referred  to  [Ref.  1].  The  specific  vector  space 


with  which  we  are  concerned  is  the  space  of  functions  defined  by 

^2[0,l]=|/:||/||=(j[1/W2^)a  <oo|. 
We  note  that  the  above  norm  is  induced  by  the  inner  product, 

(1,9)  =  f  f(x)g(x)dx. 
Jo 

That  is, 

11/11  =  ^ifJ)- 

In  particular,  we  seek  an  orthogonal  polynomial  basis,  and  derive  an 
optimization  technique  to  find  a  polynomial  pn,  of  order  n,  when  we  are  given  a 
set  of  orthogonal  polynomials,  p0,px, . . .  ,pn-i-  Recursive  application  of  a  method  for 
generating  pn  leads  to  a  complete  set  of  basis  polynomials.  The  polynomial  basis  is 
of  particular  importance,  as  the  Weierstrass  Approximation  Theorem  assures  us  that 
any  continuous  function  /,  defined  on  [0, 1],  may  be  approximated  arbitrarily  well 
with  polynomials  [Ref.  2]. 

There  are  a  number  of  existing  techniques  for  the  generation  of  or- 
thogonal polynomials.  For  example,  the  Gram-Schmidt  algorithm  may  be  applied  to 
the  sequence  {l,x,x2, . . . ,  xn, . . .}.  Another  approach  involves  solving  a  three-term 
recurrence  that  generates  the  polynomials.  We  consider  an  optimization  approach, 
in  which  we  formulate  an  optimization  problem  whose  solution  gives  us  pn.  It  is  an 
approach  that  is  suitable  for  inductive  iteration. 


2.        Example  2:    Uniform  Approximation  of  the  Expo- 
nential Function 

The  second  example  is  a  specific  problem  in  uniform  approximation. 
We  seek  the  linear  combination  of  polynomials  on  the  interval  [0,  3]  that  best  ap- 
proximates the  exponential  function  in  the  uniform  norm  sense.  The  problem,  con- 
sequently, is  to  find  the  coefficients  a,,  that  minimize  the  expression 

max  I  fit)  -  el  I, 

f£[0,3]  '      w  " 

n 

where  f(t)  =  y^  atf1. 

j=0 

We  consider  specific  cases  of  this  example.  That  is,  we  seek  the  polynomial  for  some 
fixed  degree,  n,  that  best  approximates  the  exponential  function.  Note  that  the 
uniform  approximation  problem  is  fundamentally  an  optimization  problem.  The  use 
of  the  Simplex  algorithm  to  solve  the  problem  is,  however,  unusual. 

C.      EXAMPLE    3:     THE   IMAGE    RECONSTRUCTION 
PROBLEM 

The  third  example  is  the  image  reconstruction  problem.  As  with  the  first  two 
examples,  there  are  many  existing  techniques  for  solving  this  problem.  Unlike  the 
others,  however,  this  is  an  active  area  of  modern  research,  and  the  best  methods 
of  solution  may  yet  be  unknown.  The  reader  is  referred  to  [Ref.  3]  for  a  thorough 
treatment  of  the  problem,  and  to  [Ref.  4]  and  [Ref.  5]  for  an  introduction  to  some 
recently  developed  solution  techniques. 


Suppose  a  neurosurgeon  wishes  to  rule  out  the  possibility  that  a  patient,  Fred, 
suffers  from  a  brain  tumor.  Further,  the  physician  opts  to  make  use  of  the  CAT 
(Computer  Aided  Tomography)  scan  device,  and  examine  the  inside  of  Fred's  head 
without  exploratory  surgery. 

The  CAT  scan  machine  works  by  projecting  a  finite  number  of  X-rays  of  known 
intensity  into  the  patient's  head  from  a  finite  number  of  positions.  The  intensity  of 
the  X-rays  upon  leaving  Fred's  head  is  measured.  The  intensity  of  the  emergent  X-ray 
depends  essentially  on  the  density  of  Fred's  head  over  the  locations  through  which 
the  X-ray  passes.  Having  collected  data  from  a  number  of  X-rays,  the  gathered  data 
are  processed,  forming  a  model  of  the  density  of  Fred's  head.  That  is,  the  processing 
of  the  data  results  in  the  construction  of  an  image,  and  presumably,  an  image  that 
closely  corresponds  to  the  interior  of  the  Fred's  head.  This  data  processing,  in  this 
example,  constitutes  solving  the  image  reconstruction  problem. 

1.       X-Ray  Computed  Tomography 

Understanding  the  methods  of  reconstruction  requires  that  we  know  the 
process  by  which  the  data  for  reconstruction  are  obtained.  We  begin  with  a  basic 
discussion  of  the  manner  in  which  an  X-ray  moves  through  an  object  of  homogeneous 
density,  then  derive  the  manner  in  which  it  moves  through  more  complicated  media. 

It  has  been  shown  empirically  that  the  fractional  decrease  in  beam 
intensity  of  a  narrow  beam  of  X-ray  photons  passing  through  a  homogeneous  material 


Figure  1.  An  X-Ray  Over  a  Homogeneous  Object:    I0  =  Input  Intensity.  /  =  Emer- 
gent Intensity,  p  =  Density. 


is  given  by  the  relationship  ([Ref.  6]) 


/o 


_  e-p(Ax) 


where  /o  is  the  X-ray  input  intensity,  and  /  is  the  observed  intensity  after  the  ray 
passes  a  distance  Ax  through  the  material.  See  Figure  1.  The  parameter  p  is  de- 
termined by  the  density  of  the  material.1   For  two  media,  the  fractional  decrease  is, 


predictably, 


_  _.      -(pi(Axi)+p2(Ax2)) 

I    ~  ' 

•'0 


where  Ax,  denotes  the  distance  the  X-ray  travels  through  the  i     medium 


lp  also  depends,  to  a  lesser  extent,  on  a  number  of  other  factors,  including  the  nuclear  composition 
characterized  by  the  atomic  number  Z,  a  function  of  th  ^mogeneous  material.  [Ref.  3]  pertains. 
For  the  purposes  of  this  paper,  the  effects  of  other  parai.    lers  are  assumed  to  be  nil. 


Let  us  partition  the  media  through  which  the  narrow  beam  travels  into 
n  homogeneous  segments.  Denote  the  density  over  a  single  segment  by  p(x).  The 
decrease  in  this  case  is  expressed  by 


—     =     expl-^2p{xi)Axi 
Letting  n  =>  oo,  and  Ax,  =>  0,  equation  (II.  1)  becomes 


—     =  lim       exp  l-^2p{xt)Axl 

J0  n^oo,A*-K)  y       ^ 


(II.l) 


implying 


=    exp  (  —  /  p(x)dx  )  , 


—  In  —     =      /  p(x)d: 
Iq  J 


(II.2) 


Concluding,  let  /  be  the  line  describing  the  path  of  the  X-ray,  and  the  function,  f(x,  y) 
is  the  density  of  the  media  over  /.  Let  ds  denote  a  length  over  the  line  /.  Equation 
(II. 2)  may  be  written  in  the  form 

-ln-f     =     [f(x,y)ds.  (II.3) 

Iq  Jl 

2.       The  Radon  Transform 

This  section  is  an  introduction  to  the  Radon  Transform,  and  elaborates 
its  relation  to  the  data  collected  with  the  X-ray.  We  first  define  the  transform,  then 
briefly  describe  some  of  its  properties.  The  discussion  in  this  section  pertains  to  the 
two-dimensional  case.  That  is,  we  wish  to  find  the  density  of  an  object  defined  over 
some  subset  of  7£2.  For  generalizations  into  higher  dimensions,  see  [Ref.  3]. 


We  begin  by  considering  some  density  function  /,  defined  and  bounded 
on  a  simply  connected,  compact  subset  ft  C  1Z2.  Define  L  to  be  the  set  of  all  straight 
lines  passing  through  any  portion  of  Q.  That  is,  L  =  {/  :  /  f]  ft  ^  0}.  Note  that 
the  cardinality  of  L  is  uncountably  infinite.  The  Radon  transform  is  defined  by  all 
possible  line  integrals  of  the  form: 

/  =  /  f(x,y)  ds,      j  t  J,  (II.4) 

Jij 

where  ds  is  an  increment  of  length  along  /7,  and  J  is  the  index  set  of  the  set  L. 

Consider  how  the  lines,  over  which  the  integrals  above  are  computed, 
are  determined.  Let  \i  =  [cos</>,  sin<^]T.  Then  for  a  fixed  angle  of  rotation  <f>  and 
a  distance  p  from  the  origin,  we  may  identify  each  line,  /,  by  the  set  of  vectors, 
x  =  [x,y]T,  that  satisfy  the  equation 

(x,/i)    =  x  cos  cf)  +  y  sin  (f)  =  p. 

(See  Figure  2).   Consequently,  we  may  denote  each  of  the  line  integrals  defining  the 
Radon  transform  by 

f(<t>,p)=[  /(x)dx.  (II.5) 

Again,  it  is  vital  to  note  that  the  Radon  transform  is  defined  by  the 
collection  of  all  such  line  integrals.  Consequently,  to  determine  the  Radon  transform 
fully,  we  must  know  f(4>,p)  for  all  values  of  </>  and  p.  When  we  know  the  value  of  the 
line  integrals  for  only  certain  values  of  <f>  and  p,  we  say  that  we  have  a  sample  of  the 
transform. 
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Figure  2.   The  Line,  L,  as  it  Relates  to  (p,  (f>) 

3.       The  Problem  Statement 

We  note  that  the  right  hand  sides  of  Equations  (II. 3)  and  (II. 4)  are 
identical.  We  conclude,  then,  that  if  the  X-ray  is  sufficiently  narrow,  and  we  are 
able  to  take  an  X-ray  along  all  possible  lines,  the  resultant  infinite  collection  of  data 
corresponds  to  the  Radon  transform  of  the  desired  density  function. 

The  Radon  transform  has  been  shown  to  be  one-to-one  ([Ref.  3]).  That 
is,  when  all  values  of  the  line  integrals  are  known,  one  may  determine  the  unique 
density  that  produces  the  observed  transform  data.  However,  in  most  cases  of  prac- 
tical interest,  we  are  presented  with  but  a  sample  of  the  transform  from  which  to 
reconstruct  an  image.  That  is,  we  are  able  to  collect  only  a  finite  number  of  X-rays. 
Additionally,  the  photon  beam  is  not  sufficiently  narrow  to  be  a  true  line  integral 
defining  the  transform.  In  this  case,  inverting  the  transform  is  an  ill-posed  problem. 
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If  there  exists  one  density  function  whose  sampled  Radon  transform  equals  a  given 
data  set,  then  there  exist  infinitely  many  density  functions,  /  such  that  /  =  b.  where 
b  is  the  data  obtained  from  a  transform  sample.  It  is  this  fact  that  leads  us  to 
investigate  an  optimization  approach  to  the  image  reconstruction  problem. 

D.      OPTIMIZATION 

Each  of  the  examples  can  be  formulated  as  an  optimization  problem.  Fun- 
damental to  any  optimization  problem,  and  to  the  Linear  Optimization  Problem,  in 
particular,  are  the  concepts  of  feasible  set  and  objective  function. 

1.       Feasible  Sets 

To  help  explain  a  feasible  set,  we  consider  an  example.  Suppose  we 
wish  to  model  the  production  schedule  for  a  baseball  and  softball  manufacturing 
plant.  The  company  is  required  to  make  at  least  500  baseballs  and  1000  softballs 
each  month  to  satisfy  contractual  agreements.  The  company  expects  to  procure  2,000 
pounds  of  stuffing  material,  and  3,000  square  feet  of  leather  covers.  Each  baseball 
requires  £  pound  of  stuffing,  and  |  square  feet  of  leather.  The  requirements  for  the 
softballs  are  |  pounds  and  |  square  feet  of  stuffing  and  leather  respectively.  Then 
of  all  possible  production  schedules,  we  restrict  our  attention  to  those  that  fulfill 
contractual  requirements  and  do  not  utilize  assets  which  are  not  available.  Let  6  and 
5  be  the  number  of  baseballs  and  softballs,  respectively,  produced  in  a  month.  Then 
we  require  that: 

6    >     500 

12 


s    >     1,000 

1  3 

-6  +  -s    <    2.000 
4        8       ~ 

-b  +  -s    <    2.000.  (II. 6) 

2  4       ~  y        ' 

We  have  defined  a  subset  of  all  possible  schedules  by  a  group  of  mathematical  re- 
lationships. In  this  example,  the  feasible  set  is  the  set  of  all  production  schedules 
that  satisfy  the  equations  of  (II. 6).  In  general,  we  define  the  feasible  set,  Y,  to  be  the 
collection  of  values  satisfying  the  mathematical  relationships  imposed  by  the  problem. 
2.       The  Objective  Function 

The  objective  function,  g,  defined  over  a  feasible  set,  Y,  is  the  function 
by  which  one  models  the  quality  of  a  solution.  In  the  manufacturing  schedule  example, 
we  might  logically  define  the  objective  function  to  be  profit.  Supposing  that  each 
baseball  contributes  $1  of  profit,  and  each  softball,  $.75,  we  could  write  our  objective 
function: 

g  =  b-\-  .755, 

and  we  seek  the  maximum  value  of  g  over  Y. 

Simply  stated,  an  optimization  problem  is  expressed  by  "Considering 
all  members  of  the  feasible  set,  Y,  which  member(s)  results  in  the  optimum  value  of 
the  objective  function,  glv 
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E.      LINEAR  OPTIMIZATION 

The  Linear  Optimization  Problem,  or  LOP,  is  defined  by  the  criteria  that 
the  objective  function  and  the  relationships  defining  the  feasible  set  be  linear  in  our 
decision  variables,  or  the  variables  representing  the  values  we  seek.  Then  we  may 
write  the  LOP  as  follows: 


Let  a  vector  c  =  [cj,  c2, . . . ,  cn]T  e  7£n,  a  non-empty  index  set  5,  and  for  every 
s  e  S  a  vector  a(s)  t  7£n,  and  a  real  number  b(s)  be  given.  Defining  (u,v)  as 
the  standard  inner  product,  we  seek  a  vector  y  t  7Zn,  called  the  optimal  vector, 
that  minimizes: 

(c,y) 


while  satisfying: 

(a(s),y)    >   b(s), 

for  all  s  t  S. 


We  observe  that  a  linear  maximization  problem  may  be  put  into  the  form 
above  in  the  following  way.  The  linearity  of  the  objective  function  assures  us  that  it 
is  continuous  on  Y,  and  that  the  feasible  set  is  compact.  Then  max(/)  =  min(— /), 
and  we  may  equivalently  seek  to  minimize  (— /). 

A  similar  change  may  be  made  in  the  constraints  to  reverse  inequalities  if 
necessary.  That  is,  the  problems 

Maximize:  (c,y) 

Subject  to:  (a(s),y)        <   6(5) 

for  all  s  e  S 
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and 

Minimize:  —  (c,y) 

Subject  to:        — (a(s),y)      >  —b(s) 
for  all  s  t  S 

are  identical. 

1.       The  Linear  Program 

The  case  where  the  cardinality  of  S  =  m  <  oo  defines  a  Linear  Program. 
This  special  case  of  the  LOP  is  of  particular  interest  as  it  forms  the  basis  for  finding 
solutions  to  LOPs  when  the  index  set  S  is  infinite.  Throughout  this  thesis,  the  reader 
may  assume  that  discussion  of  the  general  linear  optimization  problem  permits  the 
possibility  of  an  infinite  index  set  5,  unless  explicity  otherwise  noted. 

Now,  however,  we  examine  this  Linear  Programming  case  to  clarify  the 
concept  of  linearity.  The  problem  becomes 

minimize        (c,  y) 
subject  to:      (a(sj),y)     >  b(st) 
for  i=  l,2,...,m,     overall     y  e  Tln.  (II. 7) 

Let  aj(st)  denote  the  jth  component  of  the  vector  a(si).  We  may  write  the  problem 
as 

Minimize  cxyx  +  c2y2  H 1-  c^yn 
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Subject  to:        a^)  y1  +  a2(.Si)  y2  + h  an(si)  yn      >  b{si) 

a\(s2)  yi  +  a2(s2)  V2 -\ \- an(s2)  yn       >    b(s2) 


ai(sm)  V\  +  a2{sm)  t/2  H +  an(sm)  yn     >    b(sm) 

over  all  y  t  IT.  (II. 8) 

We  note  that  in  this  case,  we  may  define  the  feasible  set  by  the  notation 

ATy    >    b  (II.9) 

with  A  t  /JZnXm}  and  the  ith  column  of  A  is  a(s;).  The  ith  component  of  the  vector  b 
is  given  by  6(5,).  The  linearity  assumptions  can  be  expressed,  as  follows    [Ref.  7].2 

1.  Proportional :  The  objective  function  is  linear  in  the  feasible  set,  Y, 
in  the  following  sense.  Given  a  variable,  y7,  its  contribution  to  the  objective  function 
is  c3y3.  So  then  a  change  of  d  units  in  y3  results  in  a  change  in  the  objective  function 
value  of  Cjd.  Similarly,  the  constraints  are  linear  with  respect  to  the  variable  yJ? 
insofar  as  the  contribution  of  the  variable  y3  to  the  ith  constraint  is  aJ(sl)yJ.  Then 
changing  the  value  of  yj  by  d  units  changes  the  value  of  the  left-hand-side  of  the  ith 
constraint  by  aj(si)d  units. 

2.  Deterministic  :  The  components  of  the  vectors  a(s)  and  c  are  all 
determined,  as  is  each  scalar  6(5).  That  is,  if  the  components  are  derived  from  some 


2 [Ref.    7]  also  identifies  the  qualities  of  additivitiy  and  divisibility  as  requirement  of  the  linear 
optimization  problem.  These  qualities  are  deemed  to  be  inherent  in  the  qualities  defined  above. 
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stochastic  model,  their  variability  is  disregarded,  and  the  numbers  are  fixed  for  a 
given  linear  optimization  problem.  Having  defined  the  Linear  Optimization  problem, 
we  now  turn  our  attention  to  exploring  the  utility  of  solution  techniques  to  non- 
traditional  optimization  problems. 
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III.         THE  EXAMPLES,  A  DIFFERENT 

PERSPECTIVE 

A.  OVERVIEW 

This  section  addresses  some  of  the  basic  properties  of  the  sets  from  which  we 
choose  an  optimal  vector  in  our  examples.  It  is  the  structure  which  we  are  able  to 
assign  to  these  sets  that  permits  us  to  exploit  the  theories  regarding  convexity,  and 
subsequently,  the  duality  results  which  we  derive  in  subsequent  chapters.  We  then 
introduce  assumptions  that  refine  the  feasible  sets. 

B.  LINEAR  VECTOR  SPACES 

Before  proceding  to  the  specific  examples,  we  first  turn  our  attention  to  the 
matter  of  linear  vector  spaces.  A  vector  space,  L,  is  called  a  linear  vector  space  if  for 
any  vectors  x,y,z  e  L  and  any  real  scalars  a  and  a   the  following  results  hold  [Ref. 

1]: 

1.  a(x  +  y)  =  (ax  +  ay)  e  L 

2.  a(/?x)  =  (a/?)x 

3.  x  +  y  =  y  +  x 

4.  x+(y  +  z)  =  (x  +  y)  +  z 

For  each  of  the  example  problems,  the  feasible  set  is  a  subset  of  a  linear  vector 
space.  Consider  the  problem  of  finding  an  orthogonal  polynomial,  pn,  of  order  n.  It 
is  elementary  that  the  set  of  polynomials  of  order  n  form  a  linear  vector  space.  The 
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same  holds  for  the  problem  of  finding  the  polynomial  that  best  approximates  the 
exponential  function  on  [0,3].  Finally,  in  example  three,  we  have  specified  that  we 
wish  to  find  a  density  function,  /,  from  the  set  of  all  bounded,  piecewise  continuous 
functions  with  support  over  a  compact  set  Q.  The  set  of  all  such  functions  is  a  linear 
vector  space. 

Equally  important  to  our  discussion  is  the  concept  of  a  norm.  In  general,  a 
norm  on  a  linear  vector  space  L  is  defined  to  be  a  mapping,  denoted  |  |  :  L  — >  1Z+ 
satisfying  the  following  rules  [Ref.  1].  For  all  x,y  e  L,  and  a  t  7Z, 

1.  ||x||  >  0  and  ||x||  =0,<£>x  =  0 

2.  ||<*r||  =  |  a  |  ||x||, 

3.  ||x  +  y||<||x||  +  ||y||. 

Any  linear  vector  space  equipped  with  such  a  function  is  said  to  be  a  normed  linear 
vector  space.  Each  of  the  feasible  sets  of  the  examples  is  a  subset  of  a  normed  linear 
vector  space.  The  first  two  examples  are  clearly  so.  Any  norm  on  7Zn  suffice.  In  the 
third  example,  we  use  the  infinity  norm,  defined  by: 

H/IU  =  sup  |  f(u>)  I 

as  an  appropriate  norm. 
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C.      REFINING  THE  FEASIBLE  SUBSET  OF  THE  OR- 
THOGONAL BASIS  PROBLEM 

In  the  first  example,  we  are  interested  in  finding  a  polynomial,  pn,  of  order  n. 
such  that: 


(Pn,Px)      =        /     PnPidx  =  0, 

Jo 


■] 

10 

for  all  0  <  i  <  n  —  1, 


where  the  result  is  assumed  true  for  all  p,, pj,  i  ^  j.  That  is,  given  orthogonal  poly- 
nomials poiPii  •  •  ■  iPn-\i  we  seek  a  polynomial  of  order  n,  orthogonal  to  all  of  the 
polynomials  of  lower  order.  We  formulate  this  problem: 

minimize:      £"=0  /,}  pnpidx 
Subject  to:         /J  pnpi  >  0,       for  i  =  1,2, . . .  ,n  —  1.  (HL1) 

Theorem  1.  T/ie  optimal  objective  function  value  for  the  orthogonal  polyno- 
mial problem  is  zero,  and  any  optimal  vector,  pn  satisfies  the  desired  orthogo- 
nality conditions. 

Proof:  Since  we  know  triangular  families  of  orthogonal  polynomials  exist,  we  con- 
clude immediately  that  the  optimal  objective  function  value  is  bounded  above  by  zero. 
The  constraint  gives  us  zero  as  a  lower  bound.  That  any  optimal  vector  satisfies  our 
orthogonality  conditions  is  immediate  from  these  facts.  That  is,  a  zero  objective 
function  value,  in  conjunction  with  the  constraint  ensures  orthogonality.  □ 

There  are  infinitely  many  polynomials  that  satisfy  the  above  criteria.  Specif- 
ically, if  the  objective  function  evaluates  to  zero  for  some  pn,  it  clear  evaluates  to 
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zero  for  a  pn,  for  any  a  c  7Z.  Consequently,  we  add  the  additional  constraint  that  the 
polynomial  we  desire  is  the  monic  orthogonal  polynomial.  The  additional  constraint 
leads  rather  easily  to  an  n  x  n  linearly  independent  system  of  inequalities,  where  the 
unknown  element  of  7Zn  is  the  vector  whose  components  are  the  coefficients  of  the 
desired  polynomial. 

To  illustrate,  let  us  consider  the  specific  cases  of  finding  the  first  order  and  sec- 
ond order  polynomial  satisfying  (III.l).  We  input  the  zero"1  order  monic  polynomial, 
po  =  1,  to  start  the  iterative  process. 

In  the  first  order  case,  the  objective  function 

n-1      -j 

Y]  /    pn(x)  px{x)dx 
1=oJo 

is  simply 

/    pi(x)dx  =  /    (x  +  a)  dx  =  -  +  a. 
Jo  Jo  2 

The  optimization  problem  takes  on  the  form, 

minimize:      \  -\-  a 
subject  to:       |  +  a     >  0, 

from  which  we  observe  that  a  =  — |,  and  conclude  that  p\{x)  =  x  —  |.  While  the 
solution  of  this  particular  problem  is  trivial,  there  are  some  important  conceptual 
principles  working  here.  Considering  the  problem  in  terms  of  the  linear  optimization 
problem,  observe  that  the  feasible  set  is  the  set  of  all  real  numbers  a,  with  a  >  —  \. 
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As  the  function  we  seek  to  minimize  takes  on  the  form,  C  +  a,  where  C  is  a  fixed 
c   .istant,  we  clearly  wish  to  select  the  smallest  possible  value  for  a. 

Similarly,  consider  the  formulation  of  the  problem  of  finding  the  second  order 
polynomial.  We  define  the  polynomials  p0  and  pi,  as  above,  and  let  p2  =  x2-\-a-iX-\-a0. 
Computing  the  integrals,  we  find  that 


y1  f1  1      «i 

/    P0P2  =        P2  =  -  +  -r-  +  <*o, 
Jo  Jo  6        1 


and 


/    P1P2     =      I    (x--)(x2  +  a1x  +  a0)dx 


-1: 


Ql  ~  2/        v*0_T/  x~~ 


I     (SI    _\     (S!±    Sl\     Si 

4      V  3    "  6/      V  2        4/       2 
12       12 


Then  the  linear  optimization  problem 


minimize:      £?=(?  Jo  P«Pn 
subject  to:       J0  p,pn  >  0,     i  =  1, 2, . . . ,  n  —  1, 


becomes: 


minimize:      -^  +  y^ai  +  a0 


subject  to:        |  +  |qi  +  a0      >  0, 


H  +  I>         >0.  (III.2) 
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\                      Feasible  Set 

<a(s1),y>  =  -l/3V 

\ 

' 

\ 

<a(s  ),y>  =  -l/12 

a, 


Figure  3.   The  Feasible  Set  of  Example  1:  n=2 

As  we  are  currently  finding  the  feasible  set,  and  viewing  the  problem  in  terms 
of  the  general  formulation,  we  make  the  following  observations.  The  index  set  S  has 
cardinality  2.  By  rearranging  the  constraint  equations,  we  find  the  constraint  vectors 
are 


a(si)     = 
a(s2)    = 


\-  >1 

.2 

7 

\h  ' 

■jT 

with  b{s\)  —  — |,  and  b(s2)  =  —  j$.  As  the  vector  we  seek,  y  =  [o^      a0]Te  7v2,  we 
may  illustrate  the  feasible  set  as  in  Figure  3. 

We  observe  that  we  have  a  problem  of  finding  the  optimal  vector  in  7Z2  when 
seeking  the  second  order  orthogonal  polynomial.    This  property  generalizes  for  any 
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order  of  polynomial.    That  is,  if  we  seek  a  polynomial  of  order  n,  we  seek  a  vector, 
y  e  7£n,  giving  the  coefficients  for  the  optimal  monic  polynomial. 

D.      THE  FEASIBLE  SET  IN  THE  UNIFORM  APPROX- 
IMATION PROBLEM 

Consider  the  problem  of  approximating  the  exponential  function,  e*,  in  the 
interval  [0,  3]  by  a  linear  combination  of  polynomials.  We  have  specified  that  we  wish 
to  find  the  combination  that  minimizes  the  maximum  residual  over  the  interval,  and 
not  the  total  residual.  Hence,  we  are  not  solving  the  least  squares  problem,  where 
orthogonality  of  the  approximating  functions  dramatically  simplifies  the  task.  With 
the  uniform  approximation  problem,  however,  orthogonality  of  the  polynomials  is  not 
particularly  useful.  Therefore,  rather  than  using  the  orthogonal  polynomials  above, 
we  merely  specify  the  degree  of  the  approximating  polynomial.  Thus  we  seek  a  linear 
combination  of  the  polynomials 

where         />•(*)  =  t\       i  =  0,1,2, ...  n. 

Consider  the  specific  example  for  n  —  1.  We  seek  a  polynomial 

(T,y),  where  T  =  [l,t]T,  and  y  =  [a0, a^cR,2. 

Since  the  vector  T  is  fixed,  the  problem  is  equivalently  one  of  determining  the  optimal 
vector,  y  c  1Z2.  We  summarize  with  a  preliminary  statement  of  the  problem. 
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minimize:  maxte[0,3]  |  (E"=oQt^)  ~  e'  I 

over  all  y  =  [a0,  . . . ,  an]T  e  1Zn+l . 

Observe  that  the  objective  function  is  non-linear  in  the  decision  variables,  a,,  i  = 
l,...,n.  Also  observe  that  the  feasible  set  is  7Zn+1  in  its  entirety.  That  is,  any 
combination  of  real  coefficients  is  feasible,  since  there  are  currently  no  constraints. 

E.      CONVENTIONS  OF  IMAGE  RECONSTRUCTION 

For  Example  3,  we  have  specified  that  we  wish  to  find  some  function,  /,  defined 
on  the  simply  connected,  compact  set  Cl  C  7£2.  Assume  that  Q  is  a  circle  of  radius 
1.  We  also  assume  that  the  function  that  we  seek  is  piecewise  continuous  on  fi.  The 
piecewise  continuity  restriction  is  justified  by  the  physical  nature  of  the  problem  we 
are  solving.  We  call  the  space  of  such  functions  F.  Here  it  is  useful  to  define  a  basis 
for  F,  and  we  select  a  logical  basis  in  view  of  the  problem  we  wish  to  solve. 

As  we  have  stated,  the  the  formal  inverse  of  the  Radon  transform  is  well 
defined.  Our  difficulty  results  from  our  inability  to  compute  the  uncountably  infinite 
number  of  line  integrals  defining  the  Radon  transform.  This  difficulty  stems  first 
from  the  fact  that  the  region  over  which  an  X-ray  is  measured  is  not  one-dimensional. 
That  is,  the  region  over  which  the  X-ray  measures  mass  has  both  width  and  length. 
Each  X-ray  measures  the  density  of  the  medium  over  some  strip,  as  in  Figure  4. 
Additionally,  the  number  of  data  points  from  which  one  reconstructs  an  image  is 
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Source 


Detector 


Figure  4.  A  Single  Density  Measuring  Strip 


finite,  rather  than  uncountably  infinite,  as  required  for  formal  transform  inversion. 
A  more  accurate  perspective  from  which  to  view  the  data  obtained  by  the  X-rays  is 
presented  here. 

Begin  by  fixing  an  angle  <f>.  We  associate  with  each  strip  of  0,  a  label  (<f>,i). 
We  introduce  the  strip  characteristic  function,  7.  Define  the  real  valued  function  7 
defined  on  Q  by  the  rule 


7*,«M     =     < 


1,    if  a;  lies  in  strip  (<f>,  i) 
0,      otherwise. 


Then  an  integral  defining  the  sampled  Radon  transform,  for  a  fixed  angle,  <£,  and  a 
fixed  strip,  {4>,i),  becomes 


h,i  =  /  /(<*>)  7*.tM  <k>- 
Jn. 
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Figure  5.  A  Single  View:  Note  that  strips  do  not  overlap,  and  cover  Q,  completely. 


Let  us  define  a  view  to  be  the  set  of  all  strips  for  some  fixed  angle,  <f>.  We  impose  two 
restrictions.  First,  we  require  that  all  strips  of  a  view  are  non-overlapping.  Mathe- 
matically, if  (<f>,i)  and  {4>,j)  correspond  to  two  strips  of  the  same  view, 


7*,»M  70,jM  =  0,  for  all  i  ^  j, 

for  any  w  e  fi. 

Second,  we  require  that  the  strips  composing  a  view  completely  cover  the 
compact  set.  That  is,  for  any  wtfi.  and  every  angle  <j>,  there  exist  some  strip  (<f>,  i) 
such  that 

7*,i(w)  =  1. 

See  Figure  5  for  a  graphical  presentation  of  these  properties. 
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Assume  that  we  have  some  manner  in  which  to  control  the  width  of  the  strips. 
Then  we  may  select  some  number  of  strips  of  equal  width  for  each  view.  Identifying 
the  number  of  strips  for  view  <f>  as  n^,  and  the  width  of  a  strip  for  view  <p  as  6^,  we 
may  conclude  that: 

nj,  x  b^  =  1,  the  diameter  of  Cl. 

For  a  finite  number  of  views,  Nv,  the  application  of  this  convention  partition 
the  set  Cl  into  a  finite  number  of  polygons.  We  call  the  set  of  these  polygons  a 
polygonal  partition  of  Cl.  Figure  6  illustrates  the  manner  in  which  these  polygons  are 
formed.  With  each  of  the  resultant  polygons,  Sj,  we  associate  a  scalar,  area(sj),  and 
a  characteristic  function, 

(1      if  u  t  5;,  and, 
0      otherwise. 
It  is  the  set  of  these  characteristic  functions,  rftj  that  we  use  as  the  basis  for  the 
function  space,  F. 


Theorem  2.      For  any  continuous  function,  g  defined  on  Cl,  and  any  e  >  0, 
there  exist  some  polygonal  partition  on  n  polygons,  and  some  function 


f  =  Ys  ai& 
t=i 

such  that  ||/  —  gWoo  <  e. 

Note  that  we  may  write  ||/(w)  —  ^(cj)||oo  with  the  equivalent  notation, 

ll/M-0M||oo=    max   |max|  (/M-pM)^  |[. 
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Figure  6.  Polygons  Created  by  the  Views  t-^-,  {*  =  0, 1,. . .  ,4},  Each  with  4  Strips  . 

One  may  easily  verify  that  ||/(u>)  —  ^(cj)||oo  is  a  norm.  Note  that  we  may  use  the 
maximum  over  j  rather  than  the  supremum,  as  the  polygonal  partition  is  a  finite  set. 
The  properties  of  our  function  space,  F,  allow  us  to  use  the  maximum  rather  than 
the  supremum  over  each  polygon,  Sj. 

Proof:  Let  g  be  any  continuous  function  in  F,  and  let  e  >  0  be  given.  As  g  is 
continuous,  there  exists  some  8  >  0  such  that 


[x,  y)  -  (/>,  q)  Hoo     <     5  implies  that 


\\g(x,y)  -g{p,q) 


<    e. 


(III.3) 


We  use  only  two  angles,  <j>\  =  0  and  (f>2  =  \.  Let  n^  =  n^2  =  \\\-  Note  that  this 
implies  that 

n  =  rift  x  n^2  = 
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8  : 


Figure  7.  An  Arbitrary  Square  of  the  Proof  Partition 


as  in  Figure  7.     Further,  for  any  two  points,  (x,y)  and  (p,  q)  in  a  fixed  polygon, 
\\{x,y)-(p,q)\\00  <S. 

Let  /  =  JZiLi  Gifa'  We  now  consider  ||/  —  flr||oo- 

Wf-gWoo   = 


max   \  max  |  (/(u>)  -  g(u))r/>j  \ 

.7  =  1 n   {  wiSj 

f 

=      max   <  max  |  y^ajXpiip.  —  g(u)ipj 

j=l,...,n   [  wtSj        *-*> 

=      max   <  max  |  acjifrj  —  g(u>)ipj  \  > 

j=l,...,n  {  uisj  j 

=       max     <  max  I  a,  —  q(u)  I  >  . 

i=l,2 n   (.  u/eaj     '      J         »V     /   IJ 


Since  g  is  continuous  and  each  of  our  polygons  is  compact,  g  achieves  its  maximum 
and  minimum  on  each  square.  For  each  square,  Sj,  define  Mj  =  max^s  g{u),  and 
nij  =  mir^e^  g(w).  Choose 


OLj   = 


Mj  4-  m_, 
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Using  the  continuity  of  g  to  invoke  the  intermediate  value  theorem,  there  exists  some 
lo  i  Sj  such  that  g(u)  —  Qj.  Further,  we  know  that  u  c  Sj  =>  \\uj  —  uj||oo  <  S.  Therefore, 
for  any  square,  Sj, 

max  I  g(ui)  —  g(u)  I     <    £,  implying 

UltSj 

I  Qj  -  9{u)  |     <    c,  for  j  =  l,...,n. 


Therefore 


max    <^  max  |  (/(u;)  -  g(u))^j  \  \  < 

j=l,...,m   ^  utsj  ) 


D 


While  the  above  proof  uses  only  two  views,  one  may  increase  the  number  of 
views,  or  insist  on  narrower  strips  in  the  partition  of  tt.  Clearly,  such  a  refinement 
can  not  degrade  the  approximation  of  the  function  g,  but  only  maintain  or  improve 
it.  We  may,  at  worst,  maintain  the  same  constant  values  over  the  new  polygons  that 
they  were  assigned  over  the  coarser  partition. 

We  now  demonstrate  the  utility  of  defining  a  basis  for  F.  Let  k  =  YL<t>n4>- 
That  is,  let  k  denote  the  total  number  of  strips  defining  sample  transform.  For  any 
polygonal  partition  P  on  n  polygons,  the  sample  Radon  transform  with  respect  to  P, 
which  we  denote  Jp,  may  be  written  as 

fp  =  ATy 
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where  A  is  an  n  x  k  matrix,  and  y  t  7Zn.  The  matrix  A  is  given  by: 

0,  if  7,(0;)  tpi(oj)  =  0for  alio; 

At,  =  < 

area(s\),    otherwise. 
That  is,  the  AtJ  represents  the  area  of  the  ith  polygon  if  the  polygon  falls  within  strip 
j.  The  i     component  of  y  is  the  mean  density  of  the  function  /  over  polygon  i. 

For  any  fixed  polygonal  partition,  the  feasible  set  is  a  subset  of  the  infinite 
dimensional  vector  space,  F.  Each  element  of  the  subset  may  be  thought  of  as  a 
vector  in  IV1.  Without  further  restriction,  the  feasible  set  becomes  the  set  of  all 
vectors,  y  t  7Zn  such  that  ATy  =  fp.  We  exploit  many  of  the  subsequent  theorems  as 
a  result  of  the  ability  to  translate  the  problem  into  lZn. 
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IV.         CONVEXITY 

A.  OVERVIEW 

In  this  chapter,  we  investigate  the  concept  of  convexity,  both  as  it  pertains 
to  sets  and  to  functions.  The  primary  motivation  for  this  investigation  comes  from 
the  fact  that  we  may,  when  certain  convexity  conditions  are  met,  conclude  that  local 
maxima  and  minima  are  global.  Stated  differently,  we  may  eliminate  a  portion,  often 
a  large  portion,  of  our  feasible  set  from  consideration  when  attempting  to  find  the 
optimal  value  of  our  objective  function.  This  chapter  lays  the  groundwork  for  our 
investigation  into  duality,  contained  in  the  following  chapter. 

This  and  the  following  chapter  form  the  foundation  for  linear  optimization, 
and,  consequently,  the  concepts  and  results  herein  may  be  found  in  most  elementary 
texts  on  the  subject.  The  material  in  this  chapter  is  taken  primarily  from  [Ref.  8] 
and  [Ref.  9],  to  which  the  reader  is  referred  for  further  study. 

B.  CONVEX  SETS 

Let  us  return  briefly  to  the  image  reconstruction  problem.  Consider  two  ar- 
bitrary functions,  /,  g  t  F,  the  space  of  bounded,  piecewise  continuous  functions  on 
the  compact  set,  fL  Select  some  arbitrary  value  for  a  parameter,  A.  We  require  that 
A  6  [0, 1].  Consider  the  function, 

h(u>)  =  A/(w)  +  (1  -  X)g(u). 
33 


First  note  that  as  both  /  and  g  are  defined  on  ft,  so  is  h.  As  /  and  g  are  bounded  on 
a  compact  set,  M  =  max  {supn  f{u),  supn  <7(u>)}  is  well  defined.  We  know  that 

h(u)     <     AM  -f  (1  -  A)M,  implying  that 
h(u)     <    M,  for  all  u  H. 

Consequently,  the  function,  h  is  in  F.  The  important  items  to  note  here  are  that  /,  g, 
and  A  e  [0, 1]  were  each  chosen  arbitrarily.  We  conclude,  then,  for  any  two  elements 
f,g  e  F  and  for  any  A  e  [0, 1]  the  function, 

h  =  \f  +  {\-\)g     tF. 

The  above  example  proves  that  the  set  F  is  a  convex  set.  A  set  Cd,a  linear 
vector  space,  is  called  convex  if  for  any  two  elements  y,z  e  C  and  A  e  [0, 1], 

x  =  Ay  +  (l-A)z    i  C. 

Any  element  y  e  C  of  the  form  y  =  Yl\=i  My\,  with  Yl?=i  A,  =  1,  0  <  A,  <  1  is 
called  a  convex  combination  of  yi,y2,  •  •  •  yn-  This  convex  combination  is  called  strict 
if  A,  t  (0, 1)  for  all  i.  That  is,  the  convex  combination  is  strict  if  A,  ^  0  or  1,  for  all  i. 
We  now  examine  a  fundamental  characterization  of  convex  sets. 


Theorem  3.  [Ref.  8]  Let  C  be  a  convex  subset  of  L,  an  n-dimensional  linear 
vector  space.  Every  convex  combination  of  the  vectors  of  C  is  an  element  of 
C. 
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Proof:  For  n  =  1,  the  claim  is  trivial.  Assume  that  the  statement  is  true  for  r  <  ??.  —  1 
where  n  >  1.  Now  we  consider  some  convex  combination 

n  n 

y  =  Y2  x>y^  wnere  yi c  c-    H^^1*    Ai>°- 

i  =  l  i=l 

If  An  =  1,  then  we  are  done,  so  we  suppose  that  An  ^  1.  Define 

n-i  y 

A  =  y]  At-,  and  A| ■  =  — -. 

Then 

n-l 
i=l 

Note  that  sum  of  the  first  term  satisfies  the  conditions  of  the  inductive  hypothesis. 
That  is, 

n-l 

J^  a;  =  i,  and  a;  >0. 


1=1 


We  conclude  that 


'n-l 


y  =   E  Aiyi    e  c' and 


.i=l 


y  =  Ay  +  Anyn. 


Now  consider  the  expression: 


A  +  An 

n-l 


t=i 

n 

=  5> 


t=i 
=     1. 
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Then  by  the  definition  of  a  convex  set,  y  t  C.  D 

Let  A  e  7lnXTn ,  and  let  b  e  7vm.  Then  it  is  elementary  that  the  sets 

Gx     =     {x  :  ATx  =  b},  and 
G2     =     {x:ATx>b}, 

are  convex.  We  prove  the  case  of  G\. 

Proof:   Let  Xi,x2  e  Gi.  Then  Xi,x2  e  Rm,  and  ATX\  =  ATx2  =  b.    We  select  some 

value  for  A  e  [0, 1],  and  consider: 

^T(AXl  +  (l-A)x2) 
=    A^TXl  +  (l  -  A)ATx2 
=    Ab  +  (1-A)b 
=    b.  (IV.l) 

D 

One  may  show  G2  is  convex  with  an  identical  argument.  Note  that  the  set,  G2 
defines  the  feasible  set  of  the  linear  program. 

C.      HYPERPLANES,   POLYHEDRAL   SETS,   AND   EX- 
TREMA 

A  hyperplane  H  in  7ln  is  a  set  of  the  form  {y  :  (p,y)  =  k}  where  p  is  a 
nonzero  vector  in  7£n,  and  k  is  a  given  scalar.  It  is  easily  shown  that  the  hyperplane, 
H,  is  a  convex  set.    A  hyperplane  divides  7Zn  into  two  (non-disjoint)  regions,  called 
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half-spaces;  one  is  defined  by  {y  :  (p,y)  >  k}  and  the  other  by  {y  :  (p,y)  <  k,  } 
both  of  which  are  again  convex.  Note  that  the  intersection  of  a  finite  number,  m, 
of  half-spaces,  called  a  polyhedral  set,  is  also  convex,  since  the  intersection  may  be 
interpreted  as  {y  :  ATy  >  b}  where  the  i     half-space  is  define  as  the  set 

{y:(ai,y>   >  b,}. 


That  is,  A  is  an  m  x  n  matrix  whose  columns  are  the  vectors  defining  the  half  spaces. 
To  illustrate  this  point,  we  consider  a  simple  example.  Define  the  vectors 


»i  = 


,    a2  = 


1 
-1 


,  and      a3  = 


We  use  the  above  vectors  to  define  the  three  half-planes  in  7£2, 

(ai,y)>-2,    (a2,y)>-i,  and      (a3,y)  > -- 

Using  the  above  convention  for  identifying  the  matrix  A,  and  the  vector  b,  we  find 
that 


A  = 


0  1     1 

1  -1     0 


an 


d      b  = 


We  may  identify  the  intersection  of  the  half-planes  as  the  set  of  vectors,  y  in  71' 
satisfying  the  equation, 


ATy  >  b     or, 


- 

0 

-1 

y\ 

2 

1 

-1 

Vi 

> 

_i 

4 

1 

0 

1 

4 
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<a3,y>  =  l/4 
Figure  8.    The  Polyhedral  Set  of  the  Example 

The  intersection  of  these  half-planes  is  illustrated  in  Figure  8. 

We  are  interested  in  simplifying  our  optimization  problem  by  eliminating  por- 
tions of  the  feasible  set  from  consideration.  A  critical  tool  in  this  reduction  results 
from  the  notion  of  an  extreme  point.  We  here  define  an  extreme  point,  and  use  the 
concept  to  further  characterize  the  convex  sets  with  which  we  are  working  in  the 
example  problems. 

Let  C  be  a  convex  set.  We  call  y  t  C  an  extreme  point  of  the  set  C  if  it  can 
not  be  represented  as  a  strict  convex  combination  of  the  elements  of  C.  Alternately, 
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(1/4,2)                                                    / 

4  k    . « 

2 
1 

/    (9/4.  2) 

jjRk  o) 

1                                                           2                                y 

Figure  9.   Extrema  of  the  Example  Feasible  Set 

the  point  y  is  an  extreme  point  if,  and  only  if,  for  any  A  t  (0, 1),  and  for  any  x,  z  e  C, 

y      =      Ax  +  (1  —  A)z   implies  that 
y    =  x    =  z. 

Geometrically,  a  point  y  in  a  polyhedral  set  C  is  an  extreme  point  if  lies  on  some  n 
linearly  independent  defining  hyperplanes  of  C,  where  n  is  the  rank  of  matrix  AT ,  as 
formed  above.  Two  extreme  points  are  adjacent  if  the  line  segment  joining  them  is 
an  edge  of  C.  That  is,  the  line  segment  joining  them  is  formed  by  the  intersection  of 
some  n  —  1  linearly  independent  defining  hyperplanes  of  C.  See  Figure  9. 
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Theorem  4.  [Ref.  7]  Let  C  be  a  polyhedral  subset  o/7Zn.  If  C  is  bounded,  then 
C  has  at  least  n  +  1  linearly  independent  defining  hyperplanes. 

This  theorem  is  offered  without  proof.  However,  its  validity  for  the  case  of 
n  =  2  is  illustrated  in  Figure  9,  where  the  polyhedral  set  in  7Z2  has  three  independent 
defining  hyperplanes.  An  immediate  consequence  of  the  above  is  the  following: 

Theorem  5.      Let  C  be  an  arbitrary  bounded  convex  subset  o/7£n.   C  has  at 
least  n  extreme  points. 

Proof:  Suppose  that  there  are  fewer  than  n  extreme  points  of  C.  Since  any  n  linearly 
independent  hyperplanes  must  intersect  in  a  single  point  point  in  7£n,  there  are  fewer 
than  n  -f  1  linearly  independent  hyperplanes,  and  C  is  unbounded.  □ 

D.      CONVEX  FUNCTIONS 

We  now  introduce  convex  functions,  and  their  primary  characteristic  with 
which  we  are  interested.  This  introduction  is  cursory  in  nature.  For  a  more  de- 
tailed exploration  of  convexity  with  respect  to  functions,  the  reader  is  referred  to 
[Ref.  8]  and  [Ref.  10]. 

Let  C  C  7Zn  be  a  convex  set.  A  function  /,  defined  on  C,  is  said  to  be  convex 
if  for  any  elements  x,y  t  C,  and  A  e  [0, 1]: 

/(Ax  +  (1  -  A)y)  <  Af(x)  +  (1  -  A)f(y). 

If/  is  convex,  then  — /  is  said  to  be  concave.  Linear  functions  are,  thus,  both  convex 
and  concave.  Having  alluded  to  the  utility  of  convex  functions,  we  state  an  important 
result  formally. 
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Theorem  6.      [Ref.  8]  Let  f  be  a  convex  function  defined  on  a  closed  convex 
set,  C  C  TZn .   Then  a  relative  minimum  of  f  over  C  is  a  global  minimum. 


Proof:  Let  /  have  a  local  minimum  at  yi,  and  a  global  minimum  at  y2,  with  /(yi )  > 
f(y2).  Let  A  e  (0, 1)  be  given.  Because  /  is  convex 

/(Aya  +  (1  -  A)yi)    <    A/(y2)  +  (1  -  A)f(y,).  (IV.2) 

Also,  since  it  is  assumed  that  /(yi)  >  f(y2),  we  conclude 

A/(y2)  +  (1  -  A)f(yi)    <    A/(yi)  +  (l-A)f(yi) 

=  /(yi).  (iv.3) 

We  now  define  Nt(yi)  =  {y  t  Tin  :  ||y— yi||  <  e}.  That  is,  we  define  an  e  neighborhood 

about  the  point,  yi.  If 

e 
0  <  A  < -,  and, 

llyi  —  yall 

y  =  Ay2  +  (1  -  A)yj, 
then  y  e  Ne(yi).  Then 

/(y)    -    /(Ay2  +  (1-A)yi) 

<  A/(y2)  +  (1  -  A)f(yx) 

<  A/(yi)  +  (1  -  A)f(yi) 

=  /(yi), 

contradicting  IV.3,  and  the  fact  that  /  has  a  local  minimum  at  yi.  We  have  shown, 
then,  that  only  absolute  minima  are  possible.  □ 
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If  the  objective  function  is  convex  (which  it  must  be  since  we  are  considering 
only  linear  objective  functions),  we  can  be  sure  that  we  have  found  an  optimal  vector 
if  it  is  locally  optimal.  This  fact  forms  the  basis  for  the  Simplex  algorithm,  which  we 
explore  in  the  following  chapter. 

Theorem  7.  If  an  optimal  solution  to  the  Linear  Program  exists,  that  is,  if 
min  {/(y)}  exists  and  is  finite  for  some  y  in  the  feasible  set,  C  C  7£n,  then 
there  is  an  optimal  extreme  point. 

Proof:  Let  y  e  C  be  an  optimal  vector,  but  not  an  extreme  point.  Let  a  linear 
objective  function  /  defined  on  the  polyhedral  set  C  be  given.  Since  |  /  |<  oo  at 
an  optimal  vector,  one  may  clearly  add  sufficient  number  of  hyperplanes  to  bound 
the  feasible  set  if  it  is  not  already  bounded,  without  changing  the  optimal  solution. 
Assume  that  /  is  optimal  at  y.  We  consider  two  cases. 

Case  1 :  The  vector  y  does  not  lie  on  an  edge  of  C. 

We  first  recognize  that  y  can  be  written  as  a  convex  combination  of  the  extreme 
points  of  C,  since  there  are  at  least  n  linearly  independent  extreme  points.  Let 
E  =  {e  :  e  is  an  extreme  point  of  C}.  Let  E  have  cardinality  r.  Then  we  may  write 
y  =  Cj=i  ^iei-  The  linearity  of  the  objective  function,  /,  implies  /(y)  =  ZT=i  ^if(e0- 
Let  ej  be  some  extreme  point  such  that  /(e_j)  >  0,  and  let  us  decrease  the  value  of  \3 
by  5  >  0  units.  We  may  do  so  without  leaving  the  feasible  set  since  we  are  not  on  an 
edge.  Note  that  if  no  such  extreme  point  ej  exists,  then  we  may  increase  the  value  of 
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AJ?  and  the  argument  still  holds.  Call  the  new  element  of  the  feasible  set  y'.  Then 

f(y')    =    /(£Ai-ei  +  (Aj-*)ejj, 

=     /  l£A:-ei  + Ajej -Jej; 
=    /(y)-ff(e,) 

<   /(y), 

implying  that  y  is  not  the  optimal  vector,  a  contradiction.  Hence  if  y  is  a  non-extreme 
optimal  vector,  it  must  lie  on  an  edge  of  C. 

Case  2:  The  vector  y  lies  on  an  edge  of  C. 

Since  y  is  not  an  extreme  point,  but  is  on  an  edge,  it  is  on  the  line  segment 
joining  two  extreme  points,  ej  and  e2  of  C,  and  may  be  written  as  y  =  Aej  +  ( 1  —  A)e2, 
for  some  A  t  (0,1).  Parameterize  the  line  segment  between  the  points  e],e2  by  the 
equation  y(t)  =  tei  -f  (1  —  t)e2,  as  t  :  0  — >  1.  Fix  some  t  t  [0, 1],  ,  and  let  y'  = 
(1  —  t)ei  -(-  te2.  Then 

/(y')-f(y)    =  (l-0/(ei)  +  tf(e2)-Af(ei)-(l-A)f(e2) 

=  -(*-(l-A))/(ei)  +  (t-(l-A))f(e2) 

=  (<-(l-A))(-/(ei)  +  f(e2)) 

>  0,  for  all  t,  since  y  is  the  optimal  vector. 

Since  y  is  not  an  extreme  point,  it  can  be  represented  as  a  strict  convex  combi- 
nation   of    e!    and    e2.    Therefore,     we    may    choose    some    t  e  (0,  1).    such    that 
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i>  (1  -A).  Then 

t  —  (1  —  A)     >     0,  implying 
-/(ei)  +  f(e2)     >    0.  Therefore, 
/(e2)    >    /(ei). 

An  identical  argument  yields  the  result  that  /(e2)  <  f(ei).  We  conclude  that  /(ei)  = 
f(e2),  and  that  any  2  e  [0,1]  results  in  an  optimal  vector.  Choosing  t  =  0,  or  t  =  1 
places  us  at  an  extreme  point.  □ 

An  alternate  proof  may  be  found  in  [Ref.  7]. 

E.      AN  ASIDE:  THE  CONVEX  HULL 

We  desire  to  work  with  convex  subsets  of  linear  vector  spaces,  as  they  have 
useful  characteristics  when  we  attempt  to  solve  more  general  optimization  problems. 
However,  there  is  no  guarantee  that  an  arbitrary  set  is  convex.  For  such  cases,  we 
define  the  convex  hull  of  an  arbitrary  set  A  C  £,  denoted  Conv(yl),  as  the  set  of  all 
possible  convex  combinations  of  the  elements  of  A,  where  L  is  a  linear  vector  space. 
An  example  of  a  convex  hull  of  a  non-convex  set  in  7Z2  is  displayed  in  Figure  10. 

Clearly,  if  A  is  convex,  then  Conv(,4)=A  The  intuitive  notion  that  the  convex 
hull  of  a  set,  A  C  L  is  the  smallest  convex  subset  of  L  in  which  A  is  contained,  and 
conversely,  are  easily  proven  theorems  (See  [Ref.  8]). 

The  real  utility  of  the  convex  hull  stems  from  the  fact  that  any  element  of 
Conv(A)  may  be  written  as  a  convex  combination  of  the  elements  of  A.  Generating 

44 


Extensions  required  to  form  convex  hull 


Figure  10.  Forming  the  Convex  Hull 

the  convex  hull  does  not  add  any  new  extreme  points.  This  is  offered  without  proof. 
The  interested  reader  should  consult  [Ref.  8].  Consequently,  if  we  are  solving  an 
optimization  problem  with  a  linear  objective  function  on  a  non-convex  set,  ,4,  then 
solving  it  over  the  convex  hull  of  the  set  /4,  rather  than  over  the  set  A,  itself,  does 
not  change  the  solution  of  the  problem. 
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V.         DUALITY  AND  THE  SIMPLEX 

ALGORITHM 

A.      OVERVIEW 

The  concept  of  duality  makes  it  possible  for  us  to  bound  the  optimal  value  for 
the  objective  function,  /,  and  in  many  cases,  to  solve  the  LOP  more  efficiently.  As 
before,  let  c  be  a  vector  in  7Zn.  Let  S  be  an  arbitrary  index  set.  We  have  previously 
stated  that  for  every  s  t  S,  we  associcate  a  vector  a(s)  in  7£n,  and  a  scalar  b(s).  The 
general  form  of  the  linear  optimization  problem  is: 

minimize:        (c,y) 

subject  to      (a(s),y)     >   6(5),  for  all  5  e  S 

overall       y  e  Tln  (V.l) 

We  know  that  we  achieve  an  upper  bound  for  the  optimal  value  of  the  pref- 
erence function  as  soon  as  we  find  an  element  of  the  feasible  set.  However,  we  have 
no  such  simple  criteria  for  determining  a  lower  bound.  Intuitively  the  prospect  of 
finding  some  feasible  vector  is  less  daunting  than  solving  the  problem.  Using  duality 
allows  us  to  form  an  associated  optimization  problem,  find  a  feasible  vector  in  the 
associated  problem,  and  use  the  feasible  vector  to  derive  a  lower  bound  of  the  origi- 
nal problem.  The  associated  optimization  problem  is  called  the  Dual.  In  some  cases, 
we  may  bound  the  original  optimization  from  below  arbitrarily  well  using  the  dual 
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problem.  We  refer  to  the  original  linear  optimization  problem  as  the  primal.  P.  The 
primal,  P,  and  its  associated  dual,  D,  are  referred  to  as  a  Dual  Pair  . 

Define  the  value  of  a  LOP  to  be  the  optimal  objective  function  value.  We  seek 
properties  that  allow  us  to  approximate  the  solution  of  a  linear  optimization  problem 
arbitrarily  well,  and  to  determine  when  the  optimal  value  of  the  linear  optimization 
problem  and  its  corresponding  dual  are  the  same. 

This  chapter,  in  conjunction  with  the  previous  chapter,  forms  the  fundamental 
principle  underlying  the  Simplex  algorithm.  The  reader  is  again  referred  to  [Ref.  7] 
and  [Ref.  9]  for  more  detailed  descriptions  of  the  material  of  this  chapter. 

B.      WEAK  DUALITY 

We  begin  with  the  generic  linear  optimization  problem,  (V.5).  The  first  the- 
orem that  allows  us  to  bound  the  problem  from  below  is  stated  here.  Note  that  we 
allow  for  an  infinite  index  set  S. 

Theorem  8.  The  Duality  Lemma  [Ref.  9]  Let  the  finite  subset 

{S\, S2, ■ ■ ■ , sq]  C  5, 
and  the  non-negative  vector 

x  =  [xi,x2 ,xq]T 

be  such  that: 

c  =  a(si)xi  +  a(s2)x2  + 1-  a(sq)xq. 

Then  for  any  feasible  vector  y  =   [yi,y2,  •  •  •  ,yn]T   in  the  feasible  set  of  the 
optimization  problem,  P, 

b(si)xi  -\-  b(s2)x2  -r h  b(sq)xq  <  cTy. 


Proof:  Since  y  is  a  feasible  vector, 

(a(8i),y)  =  a(s,)Ty  >  b(Si),  for  i  =  1,2, . . . ,  q. 
Further,  since  xx  >  0  by  assumption, 

J2bMx*  -  H  (a(si)Ty) x" 

implying, 


t=i  t=i 


?  \T 


t=l  \i  =  l 

T 

=  c  y 


As  an  example,  consider  the  problem  of  finding  the  monic  second  order  poly- 
nomial, p2,  orthogonal  to  both,  po  =  1,  and  p\  =  x  —  |.  Recalling  from  Equation 
(III. 2),  the  primal  of  this  problem  is 

minimize    ^  +  j^a\  +  a0 
subject  to:        |  +  |ax  +  a0      >  0, 

A  +  S*i        *  °-  (v-2) 

Disregarding  the  constant  in  the  objective  function  does  not  affect  the  choice  of  an 
optimal  vector.  Consequently,  the  optimal  vector  for  (V.2)  and  the  LOP 

mininize:  (c,y) 

subject  to:       (a(s!),y)     >   6(si) 

(a(s2),y)    >  b(s2)  (V.3) 
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where 


12 


C  = 


,     a(si)  = 


,  a(s2)  = 


12 


l>M  =  -i 


an< 


b(s2)  = 


12 


are  the  same. 


Attempting  to  satisfy  the  hypothesis  of  the  duality  lemma,  we  seek  a  non- 
negative  linear  combination  of  a(si)  and  a(s2)  that  sums  to  c.  That  is,  we  seek  a 
non-negative  solution  to  the  equation: 

J'i 


1  j_ 

2  12 


1         0 


X  2 


12 

1 


Clearly,  the  only  such  vector  satisfying  the  equation  is  the  vector  x  =  [1,  1]T.  Conse- 
quently, the  optimal  value  of  the  primal  problem  of  V.3  can  be  no  better  than 

x\  b{s\)  +  x2b(s2) 

Because  the  optimal  vector  of  (V.2)  must  be  the  same  as  that  of  (V.3),  the  value  of 

(V.2)  is  bounded  below  by 

—       —    -0 
12       12  ~    ' 

as  expected. 


C.      THE  DUAL 

Having  stated  the  duality  lemma,  we  move  to  a  formal  definition  of  the  dual, 
and  similarly,  the  dual  pair.  We  begin  with  the  special  case  of  a  Linear  Program. 
The  dual  of  a  linear  program 
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minimize:  (c,y) 

Subject  to:      ATy  >  b 


is  defined  to  be 


maximize:        (b,x) 
Subject  to:    Ax  =  c 

with      xt  >  0,      fori  =  1,2,... q.  (V.4) 

Note  that  the  dual  of  an  LP  is  an  LP  itself.  To  be  feasible  above,  we  require  a 
non-negative  linear  combination  of  the  constraint  vectors  to  sum  to  the  vector  c.  The 
vector  f  becomes  the  objective  vector  in  the  dual.  These  facts  highlight  the  difficulty 
of  defining  the  dual  of  an  infinite  LOP.  Because  of  the  difficulty  of  computing  infinite 
sums  (possibly  uncountably  infinite),  we  require  a  variation  of  the  dual  for  the  infinite 


case. 


Recall  the  generic  LOP 

minimize:         (c?y) 

subject  to      (a(s),y)     >   6(5),  for  all  s  e  S 

over  all      y  t  1Zn.  (V.5) 

As  it  proves  useful  in  the  statement  of  the  dual,  we  write  (V.5)  in  the  alternate  form 


minimize  Z2r=i  ^Vr 
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subject  to:     £"=1  ar(s)yr  >  b{s) 

for  all  5  t  S.  (V.6) 

The  dual  optimization  problem,  D,  is  defined  to  be: 

Find  a  finite  subset  {si,  s2, . . . ,  sq}  C  5,  and  the  non-negative  numbers,  xi,  x2 ,  xq, 

such  that  the  expression: 

t=l 
is  maximized,  subject  to  the  constraints 

9 

5^zt-ar(s,-)     =     cv, 

for  r     =     l,2,...,n.  (V.7) 

Thaf  is,  the  dual  of  the  infinite  dimensional  LOP  is  to  find  some  optimal 
finite  subset  of  the  index  set,  and  then  solve  the  resulting  LP  dual.  In  keeping  with 
convention,  we  call  the  process  of  taking  a  finite  subset  of  an  infinite  set  discretizing.  It 
is  important  to  note  that  the  dual  is,  in  general,  a  non-linear  problem  in  2q  variables, 
since  both  the  discretization  and  values  for  the  coefficients,  x,  are  unknown.  However, 
once  we  have  chosen  a  subset  of  5,  the  problem  is  linear  in  the  unknowns  x,.  Further, 
one  might  suspect  that  if  a  sequence  of  discretizations  of  S  is  chosen  systematically, 
then  we  may  be  able  to  arrive  at  an  acceptable  approximation  of  the  solution  of 
the  associated  primal  problem,  assuming  one  exists.  That  is,  we  may  get  arbitrarily 
close  to  the  solution  of  the  dual  problem,  and  consequently,  find  an  arbitrarily  good 
approximation  of  the  solution  to  the  infinite  dimensional  primal  optimizaton  problem 
by  solving  a  sequence  of  Linear  Programs.  This  is  a  basic  premise  behind  solving 
infinite  dimensional  linear  programs  with  the  Simplex  algorithm. 
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D.      APPROXIMATING  THE  EXPONENTIAL  FUNCTION 

The  problem  of  approximating  the  exponential  function  with  an  n  degree 
polynomial  is  now  analyzed  more  closely.  Of  particular  interest  is  how  duality  results 
enable  us  to  determine  the  relative  quality  of  a  given  approximation,  and  how  they 
allow  us  to  bound  the  error  in  the  problem. 

1.       The  Primal  Problem 

Recall  that  we  stated  the  problem  of  approximating  the  exponential 
function  over  the  inverval  [0,  3]  as 

Determine  the  polynomial 


/(0  =  $>*•' 


i=0 


that  mininimizes  the  expression 


sup  |  f(t)  -  c*  I  . 

tc[0,3] 

Let  us  formulate  this  problem  in  terms  of  the  standard  linear  optimization  problem. 
We  relabel  the  index  set  T  vice  S  and  define  it  to  be  the  interval  [0, 3].  Realizing  that 
the  objective  function  above  is  a  scalar  valued  function,  as  a  first  step  we  reformulate 
the  problem  as 

minimize:  c*n+1 

subject  to:     |  Zo=o  (a^')  —  e<  I  <   an+i,      for  all  t  c  T. 

Eliminating  the  absolute  values,  we  replace  each  constraint  with  the  equivalent  pair 
of  constraints, 

n  n 

—  y]  ctjtl  +  e*    >    —  Qn+i      and,     y]  a.jtx  —  e*   >    — an+]. 

t=0  t=0 
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Rewriting,  we  arrive  at 


—  ^Q,^  +  Qn+i  >     —  et        and, 

t=0 

n 

^atf  +  an+1  >     e'. 


i=0 


Thus,  each  element  of  the  index  set  T  has  two  associated  constraint  vectors.    Let 
y  =  [q0,  ai, . . . ,  an+i]    .  We  have,  for  each  t  t  T,  a  vector 


a(t)   =     t°,  t1,  ...,tn,  1 


and  the  two  constraints 


~(a(t),y)     >     -e\  and 
(a(t),y)     >    c*. 


It  proves  useful  in  the  formulation  of  the  dual  problem  to  distinguish  the  two  con- 
straints associated  with  each  t  t  T.  As  a  notational  device  we  identify  the  vectors 

-1 
-t 
-f 


a(t+)  = 


,      and      a(t    )  = 


-tn 


It  is  important  to  note  that  the  use  of  functional  notation  for  the  vectors  a(t+)  and 
a(t~)  is  used  for  convenience  only.    No  such  functional  relationship  exists,  as  there 
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are  two  constraint  vectors  for  each  /  c  [0,3].  We  distinguish  the  vectors  by  labeling 
two  sets,  T+  and  T~ .  Note  that  T+  =  T~  =  [0,3]. 

Similarly,  for  each  t  e  T,  we  have  the  scalar,  b(t+)  =  e',  and  b(t~)  =  — e*. 
We  finally  identify  the  objective  function  vector,  c  =  [0,  0, . . . ,  1]T  e  7£n+2.  The  final 
formulation  of  the  primal  problem,  P,  is 

minimize:         c   y 

subject  to:     a(t+)Ty     >  6(*+) 

a(t")Ty     >  b(t~)  for  alU  e  T 

over  all     y  i  7ln+2.  (V.8) 

2.       The  Dual  Problem 

Having  put  the  primal  in  the  desired  form,  we  turn  our  attention  to  the 
dual.  Referring  to  the  general  form  of  the  dual  as  in  (V.7),  we  seek  the  finite  subset 
T  =  {ti,<25  •  •  •  5  tq}  C  T,  and  the  vector,  x  e  7£q,  that  maximizes  the  expression 


while  satisfying  the  constraints 

i 

^XiaT[ti)       =       Cr,  for  T  =   1 ,  .  .  .  ,  M. 


i=l 


First  make  the  substitutions  b(ti)  =  e'',  and  ar(tt)  =  i[.  As  we  have 
defined  the  set  T  to  be  T+  UT",  the  above  formulation  is  equivelent  to  the  following. 

Find  the  subsets 
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and 

f-  =  {<^2- r+}cr 

and  non-negative  scalars  x^xj, . . .  x++,  and  xf,a;J, . .  .x~  with  which  to  as- 
sociate each  element  of  the  respective  sets,  that  maximizes 

and  satisfies  the  constraints 

E*f(*?)r-E*r(<r)r    =  0,    forr  =  0,l,...,n 

t=l  1=1 

5>?  +  I>r  =  i-  (v.9) 

t=i      «=i 

The  formulation  (V.9)  may  be  written  in  the  simpler  form 

maximize:        Y11=i  €uxt 

subject  to:     Ylqt=i  x^[  —  0,      for  r  =  0, 1, . . . ,  n 

a,ii*.-i<i.  (v.io) 

where  tx  t  [0,3]  for  all  i.  The  problems  are  equivalent  in  the  respect  that  one  may 
derive  from  a  feasible  solution  of  one  a  feasible  solution  to  the  other.  The  proof  for 
this  statement  may  be  found  in  [Ref.  9]. 

3.       Qualitative  Analysis  of  Solutions 

We  begin  by  restating  the  duality  lemma  in  the  terms  of  the  uniform 
approximation  problem. 

Theorem  9.  [Ref.  9]  Let  the  finite  subset  T  C  T,  and  the  real  numbers 
X\,X2,. . .  ,xg  be  feasible  for  the  dual  problem  of  equation  (V.10).  Then  the 
following  holds  for  any  y  e  7ln+1  : 


<? 


Y2Xiet*       <       SUP 
,  =  1  tlT 
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5>r-e< 


r=0 


(V.ll) 


As  this  is  a  direct  consequence  of  the  duality  lemma,  it  is  not  proven  here,  though 
the  proof  may  be  found  in  [Ref.  9]. 

Let  us  consider  the  problem  of  approximating  the  exponential  over  T 
with  a  quadratic  polynomial.  Then  from  (V.8),  the  objective  function  vector  c  is 
equal  to  [0,0,0, 1]T.  With  each  i+and  t~  t  T  =  [0,3],  we  associate  the  vectors  and 


scalars 


a(t+)  = 


,0         i  1         ,2    „ 
t+       t+       t+       1 


andb(t+)    =   et+, 


and 


a(r)  = 


U  1  A 

-t"    -t"    -t"  1 


and  b(t    )   =  e  l 


respectively.  The  dual  problem,  from  equation  (V.10),  is  to  find  the  set  {t\,  ti, . . .  ,tq} 
T  C  [0,3]  and  associated  non-negative  scalars  that  maximize 


while  satisfying 


I>.e'"' 


i=l 


y]  Xjtj     =    0,  for  r  =  0,1,2,  and 
t=i 

ENI  <  i. 


(V.12) 


t=i 
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Let  us  arbitrarily  choose  the  subset  T  to  be  {0, 1,2,3}.  Hoping  to  apply  the  restated 
duality,  Theorem  9,  we  first  require  a  solution  to  the  equation: 


1111 
0  12  3 
0     14    9 


■ 

X\ 

0 

x2 

= 

0 

X3 
£4 

0 

(V.13; 


Every  such  vector  is  of  the  form  x  =  [—a,  3a,  —3a,  a]T,  where  a  is  an  arbitrary 
real  numbe:.     Scaling  in  order  to  satisfy  the  constraint,  ]C*-i       xi    l<    1,  we  let 

iT 


X  = 


1      3      _3      1 

81     8'         8'     8 


.  The  hypothesis  of  the  Theorem  9  satisfied,  we  conclude  that 
the  best  quadratic  approximation  to  the  exponential  function  over  T  =  [0,  3]  in  the 
uniform  norm  sense,  differs  from  ef  by  at  least: 


5>v  = 


_Ie0,    3    1  _3    2    ,     1    3 
8  '8  8        ~  8 


.6340. 


«=i 


E.      STRONG  DUALITY 

Consider  the  three  different  possibilities  we  may  encounter  in  the  solution  of 
the  Linear  Optimizaton  Problem.  Referring  to  the  optimal  objective  function  value 
of  the  minimization  problem  as  V{P),  and  to  the  optimal  value  of  the  dual  as  V(.D), 
we  list  the  possible  conditions,  or  states,  of  the  problem  as  follows  [Ref.  9]: 


Inconsistent:    (IC)     The  feasible  set  is  empty,  so  that  no 
solution  is  possible. 

Bounded:         (B)      There  exist  at  least  one  feasible  vector,  and 
among  such  feasible  vectors,  at  least  one  is  optimal. 
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Unbounded:     (UB)   There  are  feasible  vectors  such  that 

the  objective  function  may  be  made  arbitrarily  small. 

A  duality  gap  is  said  to  occur  when  V{P)  ^  V(D),  that  is,  when  the  optimal 
values  of  the  dual  pair  are  not  the  same.  We  hope  to  find  general  conditions  that 
preclude  the  existence  of  a  duality  gap.  Theorems  that  allow  us  to  disregard  the 
possibility  of  a  duality  gap  are  called  strong  duality  theorems. 

1.       The  Dual  and  Convexity 

We  briefly  characterize  the  dual  problem  as  it  relates  to  our  discussion 
of  set  convexity.  Before  continuing,  we  require  the  definition  of  the  Convex  Cone.  Let 
C  be  a  convex  subset  of  7ln.  The  convex  cone  of  C,  denoted  x{C),  is  defined  to  be 
the  set  of  all  vectors  y  e  7£n,  such  that  y  =  Ax,  where  A  >  0,  and  x  t  C. 

In  Chapter  IV  we  constructed  an  example  of  a  polyhedral  set  using  the 
vectors 


a(si)  = 


,    a(s2)  = 


-1 

1 


,      and  (s3)  = 


The  resultant  polyhedral  set  is  illustrated  in  Figure  8.  The  darkened  region  of  Figure 
11  illustrates  the  addition  to  the  set,  that  together  with  the  original  polyhedral  set, 
forms  the  convex  cone.  The  darkened  portions  of  the  Figure  extend  to  infinity. 

Consider  the  specific  case  of  the  convex  cone  of  the  constraints  of  the 
linear  optimization  problem.  We  have  expressed  the  constraints  by  (a(s),y)  >  b(s), 
for  all  s  e  S.  Define  As  =  {a(s)  :  s  e  S}  C  7£n.  We  know  that  As  is  convex  from 
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<a3,  x>=  1/4 


Figure  11.  Formation  of  the  Convex  Hull. 

equation  (IV.  1).     We  refer  to  the  convex  cone  of  As  as  the  moment  cone  of  the 
optimization  problem,  P,  and  denote  x{As)  by  Mn. 

Having  defined  the  moment  cone,  we  arrive  a  fundamental  characteri- 
zation of  the  dual  problem,  D. 


Theorem  10.  [Ref.  9]  The  dual  problem,  Z),  is  feasible  (i.e.  the  feasible  set  is 
not  empty)  if  and  only  if  c  t  Mn. 


The  proof  may  be  found  in  [Ref.  9].  The  result  follows  directly  from 
the  definition  of  the  dual.  An  alternate  interpretation  of  this  result  is  as  follows.  The 
dual  problem  is  feasible  if  and  only  if  we  may  express  the  vector  c  as  a  non-negative 
combination  of  the  constraint  vectors  of  the  linear  optimization  problem,  P. 

The  following  is  a  generalization  of  the  theorem  that  allows  us  to  express 
every  element  of  the  convex  set,  A,  as  a  convex  combination  of  the  extreme  points. 
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The  theorem  proves  vital  in  the  discussion  of  the  Simplex  algorithm,  as  it  allows  us 
to  bound  the  required  number  of  elements,  sq  i  S,  in  the  discretization  of  our  index 
set  when  forming  the  dual. 

Theorem  11.     The  Reduction  Theorem  [Ref.  9]  Let  the  vector  z  t  7Zn  be 
a  non-negative  linear  combination  of  the  vectors,  Zi,Z2, . . .  ,  zq.  That  is, 

q 

Z  =  X^XiZi' 
i=l 

with  xx  >  0  for  all  i.  Then  we  may  also  write: 

q 
z  =  ^2x'i  z''  with  X;  >  0, 

i=l 

where  at  most  n  of  the  numbers  x\  are  nonzero.   Moreover,  the  set  of  vectors 
{z;}  corresponding  to  the  nonzero  scalars  x\  are  linearly  independent. 

Proof:  We  first  note  that  if  Zi,Z2,...,zq  are  already  linearly  independent,  then 
q  <  n,  and  the  initial  representation  of  z  already  satisfies  the  theorem.  Assume, 
then,  that  q  >  n,  and,  consequently,  that  the  vectors,  Zi,Z2,...,zq  are  not  linearly 
independent.  Then  we  know  that  we  may  write 


t=i 

where  at  least  one  at  j£  0.  For  any  r  :  ar  ^  0,  we  have: 


zr=-E^i-  (V.14) 


Substituting  into  the  equation  of  our  hypothesis,  we  have: 


z  =  Yl  (Xi  _Xr~)z' 


Qr 
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We  have,  then,  a  representation  of  z  by  a  linear  combination  of  q  —  1  of  the  vectors, 
Zj.  We  must  show,  then,  that  the  expression  (xi  —  i^)  may  be  made  non-negative, 
for  i  =  1,  2, . . . ,  r  —  1,  r  +  1, . . . ,  q.  Select  some  ar  >  0.  We  can  clearly  do  so,  as  if 
all  q,  are  negative,  we  may  multiply  by  —1  and  still  have  the  desired  result  that 
X3t=i  —atz,  =  0.  Then  in  equation  (V.14),  if  a,  <  0,  we  may  conclude  that 

x,  —  xr —  >  U, 

since  x,,  j      and  ar  are  each  nonnegative. 

We  now  consider  the  case  that  a,  >  0.  Then  we  must  show  that  -^  >  — . 

*     —  Q,      —     Qr 

We  may  accomplish  this  quite  simply,  by  selecting  the  r  that  minimizes  the  expression, 
2*-  over  all  ar  >  0.  We  have  expressed  z  as  a  non-negative  linear  combination  of  q  —  1 
of  the  vectors,  Zi,z2, . . .  ,zq,  and  may  continue  inductively  until  we  have  the  desired 
result.  □ 

The  reduction  theorem  yields  this  immediate  result.  Let  S  =  {si, . . .  ,s9}  C 
S,  and  the  set  of  non-negative  numbers  {xi, . . .  ,xg}  be  feasible  for  the  dual  problem, 
D.  That  is: 

^xtar(st)     =     cv, 
t=i 

for  r  =  1,2, ...  ,n.  Then  there  is  a  subset,  S'  =  {stl  ,...,$,-„}  and  a  set  of  non-negative 
numbers,  {x't  ,...,x't  }  that  is  also  feasible  for  D.  Note  that  we  have  not  included 
the  objective  function  of  the  dual  in  our  reduction  above.  It  is  not  necessarily  true, 
then,  that  we  need  only  to  consider  discretizations  of  S  with  cardinality  n.  That  is, 
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let  us  reduce  the  non-negative  linear  combination 


2J  xt-a(sj),  where  q  >  n 


1=1 


to  the  combination 


XX-a(si), 


t=i 


where  no  more  than  n  of  the  scalars,  x'  are  non-zero.  Then  it  may  be  that 


^2xib(si)  ±  Y^x'^&i). 


1=1 


1=1 


Consequently,  we  include  the  optimal  objective  function  value  in  the  set  of  equations 
for  reduction.  This  convention  requires  that  we  define  a  new  moment  cone,  which  we 


call,  Mn+\. 


Mn+i  =  x(A') 


where  A'  is  formed  by  the  vectors, 


a'(s)  = 


b(st) 
ai(si) 


e7ln+1,  i=  1,2,. ..,n. 


The  dual,  then,  may  be  stated 


maximize:    cq 


subject  to:     c     cMn+i,  where  c  =  [c0,Ci, ...  ,cn] 


This  formulation  is  useful  in  discussion  of  strong  duality  results. 
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2.        Solvability  Conditions 

We  move  from  the  infinite  case  to  the  case  of  a  finite  index  set.  The 
following  results  are  presented,  without  formal  proof,  though  they  may  be  found  in 
[Ref.  9]  or  [Ref.  10].  These  theorems  enable  us  to  determine  when  the  dual  problem 
has  a  solution.  That  is,  we  seek  to  determine  when  there  exists  at  least  one  vector  of 
our  feasible  set  that  minimizes  our  objective  function.  Note  the  distinction  between 
solvability  and  boundedness  as  defined  in  the  state  section  above.  That  is,  we  may 
have  feasible  vectors,  but  no  optimal  vector  in  our  feasible  set.  The  discussion  in 
this  section  pertains  to  the  finite  case  of  the  linear  optimization  problem.  Readers 
interested  in  an  examination  of  some  criteria  for  the  convergence  of  the  LOP  in  the 
case  of  an  infinite  index  set  are  referred  to  [Ref.  11]. 

Theorem  12.  [Ref.  9]  Let  the  linear  optimization  problem,  P,  be  such  that 
Mn+i  is  closed,  and  the  dual  problem,  D,  is  bounded.    Then  D  has  a  solution. 

The  proof  of  this  theorem  is  straightforward.  Recognize  that  the  objective  function 
of  the  dual  is  /  :  7Zn+1  — >  TZ  by  f{zQ,Zi, . . .  ,zn)  =  z0.  Then  /  is  clearly  continuous, 
on  a  compact  set,  and  we  conclude  the  result. 

Theorem  13.  [Ref.  9]  Any  convex  cone  P  defined  by  a  finite  number  of  vectors 
in  lZn  is  closed,  in  that  any  convergent  sequence  of  vectors  in  P  converges  to 
a  vector  in  P. 

Coupling  these  observations,  we  conclude  that  any  finite  dual  pair, 
(P,  D),  with  both  P  and  D  consistent,  is  solveable.  That  is,  both  the  primal  and 
dual  have  solutions. 
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H(y,0) 

=  {  x:  <x,y>=0} 


Figure  12.   The  Separating  Hyperplane  H(y,v)  of  the  Set,  M,  at  the  Point  z. 
3.       Separating  Hyperplanes 

We  now  address  the  final  tool  that  we  use  to  eliminate  a  duality  gap 
in  the  linear  program.  Let  H(y,v)  =  {x  c  1ZP  :  yTx  =  u}.  Then  the  hyperplane, 
H(y,is),  is  said  to  separate  z,  a  vector  not  in  M,  from  the  convex  set,  M,  if 

y  x  <v<  y'z, 

for  all  x  t  M.  Figure  12  illustrates  one  seperating  hyperplane  between  the  point  z 
and  the  set  M  which  is  contained  in  7Z2.  Let  Zo  be  the  vector  in  M  closest  to  z  in 
the  Euclidean  norm  sense.  Let  y  =  z  —  Zo,  and  let  v  —  0.  Then  H(y,v)  is  the  line 
orthogonal  to  y  at  the  point  z0. 


Theorem  14.  The  Separation  Theorem  [Ref.  9]  Define  ||x||  to  be  standard 
Euclidean  2-norm.  Let  M  C  7lp  be  a  non-empty,  closed  convex  set,  and  let  z 
not  be  in  M.  Further,  let  Zq  be  the  unique  vector  in  M  such  that  ||z  —  Zq||     < 
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||z  —  x||  for  all  x  c  M.1  Finally,  let  y  =  z  —  z0,  and  1/  =  (z  —  z0)Tz0.   TTjo)  £/ie 
hyperplane,  H{y,u)  separates  z  /rom  M. 


Proof:  Let  xtM,  and  fix  0  <  /j  <  1.  Then 

(1  -  /i)z0  +  [ix  =  z0  +  /x(x  -  z0)  c  M, 

as  M  is  a  convex  set.  Further, 

||z-z0||2     <     ||z  -  (z0  +  ^(x- z0))||2 

=     ||z  -  Zo||2  -  2/*(z  -  z0)T(x  -  z0)  +  /i2||x  -  z0||2, 


which  implies  that 


(z-z0)T(x-z0)    <    -^||x-z0||2. 


Let  fj,  — >  0.  Then 

(z  —  z0)Tx     <     v,  for  any  x  e  M. 

Then  by  the  definition  of  a  separating  hyperplane,  we  have  only  to  show  that  u  <  yTz. 
Since  z  is  not  in  M, 

0<  ||z-z0||2  =  (z-z0)T(z-z0) 

T  T  T 

=  y1z-y1z0  =  y1z-K 

D 

The  separating  hyperplane  defined  above  is  a  necessary  tool  in  the  elim- 
ination of  duality  eaps  in  the  finite  linear  optimization  problem. 


'That  such  a  unique  vector  exists  is  proven  in  [Ref.  9]. 
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4.        The  Strong  Duality  Theorem 

We  close  this  section  with  a  statement  and  proof  of  a  fundamental 
theorem  of  linear  optimization,  which  states  sufficient  conditions  for  the  absence  of  a 
duality  gap  in  the  dual  pair,  (P,  D). 

Theorem  15.    [Ref.   8]  Let  the  dual  pair,  (P,D)  satisfy  the  following  assum- 
tions. 

1.  The  dual  problem  is  consistent  and  has  a  finite  value  V(D). 

2.  The  moment  cone,  Mn+\  is  closed. 

Then  (P)  is  consistent,  and  V(P)  =  V(D).  That  is,  no  duality  gap  occurs. 

Proof:  Let  let  the  vector,  c  =  [co,  Ci, . . . ,  cn]T  t  Mn+i,  be  an  optimal  solution  of 
the  dual  problem.  Then,  for  any  e  >  0,  the  vector,  c'  =  [co  +  £,  c1? . . .  ,cn]T  is 
not  in  Mn+\.  As  we  are  assuming  that  Mn+1  is  closed,  we  conclude  that  there  is 
a  hyperplane  separating  the  vector  c'  from  Mn+i .  Consequently,  there  exists  some 
vector  y  =  [yo,yi,  •  •  •  ,yn]T  <■  7£n+1,  with  y  /  0,  such  that 

n  n 

Y^xrVr    <   0   <   y0(c0-\-e)-r^2cTyr, 

r=0  r=l 

for  x  =  [xo,  Xi, . . . ,  xn]T  t  Mn+1.  Let  x  =  c.  Then  y0  e  >  0,  implying  y0  >  0.  Now  let, 

x  =  [b(s),ai(s),...,an(s)]T  e  A;   C   Mn+1. 

where  s  c  S,  and  ar(s)  is  the  rth  component  of  the  constraint  vector  associated  with 
s.  Then 

^i        \  yo  J 
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which  implies, 

T 


i 


y  = 


-yi 


1  •  •  •  1 


yo  yo 

is  feasible  for  the  primal,  P.  Further, 


i  IT 


0  <  y0(c0  +  c)  +  ^2cryr. 


r  =  l 


=>    y^Cr        L         <    Co  +  £• 

Applying  the  duality  lemma,  we  conclude  that 


V(P)   <  'Ecd  <  co  +  £  =   V(D)  +  e  <  V(P)  +  £, 

r=l 

implying 

V{P)-e<  V{D)  <  V{P), 

for  any  e.  D 

Of  a  final  note,  if  the  index  set  of  our  constraints  is  finite,  then  we  may 
conclude  immediately  that  no  duality  gap  exists  in  the  dual  pair,  (P,  D).  This  follows 
directly  from  the  above  theorem  in  conjunction  with  Theorems  12  and  13. 

F.       THE  SIMPLEX  ALGORITHM 

We  present  a  very  brief  introduction  to  the  Simplex  algorithm,  and  use  it  to 
solve  a  simple  LOP.  This  section  is  not  intended  to  illustrate  the  implementation 
of  the  algorithm  in  any  specific  form.  Rather,  this  section  attempts  to  explain  the 
algorithm  as  it  exploits  the  results  of  the  duality  concepts  above.  The  problem  is 
assumed  to  be  infinite-dimensional  in  this  presentation. 

67 


We  begin  with  a  problem,  P,  of  the  form: 

Minimize:         ]Cr=i  Cr^ 

subject  to:       5Zr=i  ar{s)yr     >    ^(s)i  f°r  aU  s  e  S. 
Then  we  write  the  dual,  D: 

Maximize:       Zl,n=i  6(3t-)xt- 

subject  to:      Y^=iar{sz)xl    =  cv,r=  1,2, ...,n 
6j  e  5,  x,  >  0. 


T 


Choose  some  subset,  cr  =  {si,s2,  ■  ■  ■ ,  5n}  C  5,  and  a  vector  x  =  [xi,X2, . . .  ,xn] 
that  is  feasible  for  the  dual.  The  methods  to  arrive  at  an  initial  feasible  vector, 
provided  one  exists,  may  be  found  in  any  Linear  Program  text.  In  particular,  the 
reader  is  referred  to  [Ref.  7].  We  derive  a  vector  y  from  our  choice  of  <j,  which  is 
associated  with  the  primal  problem.  As  a  matter  of  convenience,  we  abbreviate  this 
set  of  values,  {cr, x,y}.  We  require  that  the  vectors,  a(s;)  be  linearly  independent. 
That  we  may  always  find  a  set  of  linearly  independent  vectors  is  assumed  in  this 
presentation. 

Forming  our  matrix  A  as  before,  we  know  that  the  linear  independence  of  the 
vectors  ensures  that  there  is  a  unique  vector,  x  satisfying: 

Ax  =  c, 

since  we  are  feasible  in  the  dual.     Define  the  discretized  primal  to  be  the  linear 
program  that  results  in  considering  only  the  finite  subset  of  the  index  set  S.  Let 
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A(s\,  s2,  ■  ■  ■ ,  5n)  =  [a(s] ),  a(s2 ),...,  a'  sn)],  with  b( Si, ...,  sn)  defined  in  the  same  'ban- 
ner. From  the  discretization,  cr,  we  .  ,ok  for  a  vector,  y,  that  is  feasible  for  tl  lis- 
cretized  primal,  P.  We  note  that  one  such  vector,  y  solves  the  equation: 

AT(sus2,...,sn)y  =  b(si,s2,...,sn). 


Then 


T\-li 


y  =  (A1)"1b 


The  set  of  values  of  a  and  the  vector  y  that  is  formed  in  the  manner  above  is  called 
a  basic  solution  of  the  LOP.  The  steps  of  the  algorithm,  to  this  point  are: 

1.  Select  a  subset,  <r  C  5,  such  that  the  vectors,  a(si ),  a(s2), . . . ,  a(sn) 
are  linearly  independent. 

2.  Compute  the  unique  non-negative  solution  to  the  equation.  Ax  -  :*. 

3.  Compute  the  solution  to  the  system,  ATy  =  b, 
for  the  discretized  primal. 


Return  to  the  problem  of  approximating  the  exponential  with  a  quadratic  poly- 
nomial er  the  interval  [0,3].  We  have  formulated  the  problem  with  the  constraint 
vectors  of  the  index  sets,  a(t+),  and  a(t~),  given  by 


a(t+)  = 


1 


,     and      a(t   )  = 


-1 
-t 
-t2 

1 


Additionally,  the  constraint  scalars  were  defined  to  be  b(t+)  =  e\  and  6(2    )  =  —  e', 
and  the  objective  function  vector  was  given  by  c  =  [0,0,0, 1]T.  The  problem  is 
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minimize:         c  y 


subject  to:     a(t+)Ty     >  b(t+) 


a(t-)jy     >  b(t~)  for  all*  e  T 


over  all      y  t7ZA. 

Step  One:  Arbitrarily  select  a  to  be  composed  of  the  union  of  the  sets  V\  = 
{0,2}  CT"  and  <72  =  {1,3}  C  T+. 

Step  Two:  Compute  the  solution  of  the  system 


■ 

1   1 

-1   1 

Xi 

0 

0    1 

-2    3 

2*2 



0 

0    1 

-4    9 

%3 

0 

1   1 

1     1 

X4 

1 

(V.15) 


The  solution  of  this  system  is  given  by 


x     = 


13  3  1 

8  8  8  8 


Step  Three:  Compute  the  solution  of 


-10        0     1 


1111 


-1     -2    -4     1 


9    1 


■ 

y\ 

-e° 

V2 

— 

e1 

ys 

-e2 

k 

e3 

(V.16) 
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Figure  13.  A  First  Approximation  of  the  Exponential  Function. 

The  vector  y  =  [1.6342    —2.2946  2.7445  .6342]     is  the  unique  solution  of  this  system. 
That  is,  y  is  feasible  for  the  discretized  primal.  The  first  approximation  is  given  by 


p2(x)    =    1.6342  -2.2946X  +  2. 7445x: 


(V.17) 


The  graph  of  the  exponential  versus  the  approximation  is  given  in  Figure  13. 

We  here  introduce  a  lemma  that  offers  us  a  termination  criteria  for  the  algo- 
rithm. 


Theorem  16.  The  Complementary  Slackness  Theorem  [Ref.  9]  Let 
the  set,  {(j,  x,y}  be  as  above.  If  the  vector  y  is  feasible  for  the  non- discretized 
primal  P,  and  the  following  holds: 

x«  (  Har(s«)z/r  -M5«))  =  0,  for  r  =  1,2,..., n. 

Then  we  may  conclude,  that  if  the  vector,  y,  as  determined  in  step  3,  is  feasible 
for  the  primal,  P,  we  have  found  the  optimal  vector  in  our  problem. 
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Figure  14.  Absolute  Error  in  the  First  Exponential  Approximation. 

In  the  current  approximation  problem,  we  find  that  the  current  solution  does 
not  satisfy  this  criteria.  We  observe  the  graph  of  the  absolute  difference  between  the 
functions  and  find  that  the  error  exceeds  .6342  over  the  latter  portion  of  the  interval. 
See  Figure  14. 

The  remainder  of  the  algorithm  is  a  sequence  of  exchange  steps  that  replace 
existing  elements  of  the  set,  cr,  with  elements  that  improve  the  value  of  the  dual 
problem,  Z),  and  consequently,  improve  the  bound  of  V(P).  The  method  of  selecting 
new  elements  to  the  set,  cr,  may  change  with  implementation,  but  it  should  be  noted 
that  exactly  one  element  of  the  set  a  is  replaced  at  a  given  step,  in  any  implementation. 
Recalling  from  our  discussion  of  extreme  points  of  our  feasible  set,  that  strategy 
ensures  that  the  algorithm  looks  to  adjacent  extreme  points  for  optimality. 

We  conduct  one  such  exchange.    Note  that  the  error  is  most  severe  at  the 
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Figure  15.  Absolute  Error  in  the  Second  Exponential  Approximation. 

point  t  =  3.  Then  it  is  logical  to  seek  a  better  solution  at  that  point.    Then  we  let 
o~\  —  {0,3}.  The  new  system  of  equations  requiring  a  solution  in  step  3  is 

-10        0     1 


1111 


1     -3    -9     1 


13       9     1 


■ 

V\ 
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Vi 

— 

e1 

V3 

-e3 

k 

e3 

(V.18) 


The  solution  of  the  above  system  is  given  by 


o,o,  J,  \ 


We  find  that  the  error  is  decreased.  The  absolute  error  is  given  in  Figure  15. 

We  note  that  the  solution  is  not  feasible  for  the  entire  interval,  since  there 


exist  points  where  the  error  exceeds  .5.  Thus,  we  would  look  to  adjacent  extreme 


point  solutions  and  repeat  the  process  until  we  arrive  at  a  discretized  solution  that  is 
feasible  throughout  the  interval. 
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VI.         RECONSTRUCTION  FORMULATION 

AND  SOLUTION 

A.  OVERVIEW 

Having  laid  the  complete  foundation,  we  formulate  the  image  reconstruction 
optimization  problem.  The  first  portion  of  this  chapter  addresses  the  conceptual 
aspects  of  the  problem,  while  in  the  latter  portion  we  use  the  Simplex  algorithm  to 
solve  a  simple  reconstruction  problem.  We  conclude  the  chapter  with  a  brief  discussion 
of  the  merits  and  drawbacks  of  a  Linear  Programming  approach  to  the  reconstructs 
problem. 

B.  TARGET  FUNCTIONS  AND  NORMS 

The  problem  we  wish  to  solve  is  to  find  the  density  function,  /,  that  produces 
the  observed  sampled  Radon  transform.  As  the  problem  is  ill-posed,  we  must  define 
some  preference  function  by  which  to  compare  the  quality  of  the  infinitely  many 
density  functions  that  satisfy  the  above  requirement.  We  do  so  by  specifying  some 
function,  g,  defined  over  Q,  which  is  assumed  to  represent  the  most  likely  density 
of  the  image.  That  is,  of  all  density  functions  that  produce  the  observed  transform 
data,  we  seek  that  which  is  most  like  what  we  expect  to  find.  How  we  determine  the 
function  g  is  not  a  matter  of  discussion  here.  We  only  assume  that  we  know  some 
such  function. 
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The  problem  of  how  to  determine  the  best  solution  becomes  one  of  finding 
the  density  function  that  produces  the  observed  transform  that  is  "closest1"  to  g  in 
some  sense.  We  choose  the  infinity  norm,  or  max  norm,  to  measure  closeness.  Let 
P  be  a  polygonal  partition  of  the  compact  set  17  C  7£2,  consisting  of  the  n  polygons, 
Si,  52, . . . ,  sn.  Recall  that  the  function  i^j{ui)  is  defined  to  be  the  characteristic  function 
of  polygon  j  in  P.  Imposing  the  restriction  that  the  optimum  density  be  constant  over 
each  polygon,  the  density  takes  on  the  form, 

n 
3=1 

We  seek  a  density,  /  defined  over  £7  that  minimizes  the  following: 

ll/M-sMlloo     =       max    |max|ai-flr(a;)|[.  (VI. 1) 

.7  =  1,2,.. .,n   ^  ujcSj  J 

We  also  choose  some  i  >  0  and  insist  that 

fp  —  t>     <     £3  and 
/    >    0, 

where  the  vector  inequality  is  componentwise.  Recall  that  fp  is  defined  to  be  the 
sampled  transform  of  the  density  /  for  partition,  P.  The  vector  b  is  the  observed 
sample  Radon  transform.  The  non-negativity  constraint  stems  from  the  physical 
nature  of  the  problem.  That  is,  we  do  not  accept  solutions  that  attribute  negative 
density  to  physical  objects. 

Before  continuing,  let  us  consider  the  objective  function  of  equation  (VI. 1). 
Recall  that  our  attention  is  fixed  on  density  functions  defined  on  fi,  a  compact  subset 
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of  7v2.  Let  us  first  fix  our  attention  on  some  polygon,  Sj,  in  the  polygonal  partition  P. 
Let  Mj  denote  the  largest  absolute  difference  between  our  target  function,  g  and  the 
scalar,  ctj,  that  we  associate  with  the  polygon,  Sj.  That  is,  f(u)  —  Qj,  for  all  cj  e  sr 
The  term 

max  |  atj  -  g{u)  \ 

UltSj 

is  well  defined,  as  both  functions  are  piecewise  continuous  over  the  compact  set,  Sj. 
The  objective  function  is  defined  to  be  the  largest  of  the  M3  values  over  all  polygons. 
As  the  problem  is  not  linear,  we  write  an  equivalent  formulation: 

minimize:  k 

subject  to:     \\fp  —  b||oo  <    i 

ctj  -\-  k  >     g(u)rjjj(uj),       for  all  u  e  fi, 

—atj  +  k  >     —g(u)xl)j(u>),       for  all  to  e  $7, 

a,  >  0,       for  all  j.  (VI.2) 

Suppose  the  target  function,  g,  is  chosen  to  be  continuous  over  ft,  and  further 
suppose  that  gp  =  b,  where  gp  is  the  sample  transform  of  the  target  density,  g. 
If  the  method  is  to  prove  worthwhile,  we  expect  that  the  test  density  function  is 
optimal.  That  is,  if  the  test  and  target  densities  are  the  same,  we  can  expect  to  find 
an  arbitrarily  good  approximation  of  the  test  density.  We  state  the  above  formally. 


Theorem  17.       Let  g  be  a  non-negative,  continuous  function  defined  on  the 
set  ft.  Additionally,  let  values  for  e  >  0  and  i  >  0  be  given.    Then  there  exists 


1 1 


some  partition,  P,  of  n  polygons,  and  an  associated  function,  f  =  ]C?=i  ^j^jt 
so  that  the  optimum  value  of  the  linear  optimization  problem: 


minimize:  k 

subject  to:      \\fp  -  gp\\oo  <     S  (VI. 3) 

CLj  +  k  >     g{uj)ipj,      for  all  u  t  fi,  (VI. 4) 

—ctj  +  k  >     —g(uj)\l)j,      for  all  u  t  Cl,  (VI. 5) 

Qj  >  0,       for  all  j  (VI. 6) 

is  less  than  e. 

Proof:  We  show  that  we  may  find  a  feasible  vector  for  any  value  of  fc,  and  con- 
sequently, for  k  <  e.  The  proof  depends  on  the  continuity  of  the  sample  Radon 
transform.  That  is,  let  g  be  any  continuous  function  defined  on  ft,  and  let  i  >  0  be 
given.  Then  we  wish  to  show  that  there  exists  some  Si  >  0  and  some  partition  Ps1 
such  that  the  following  property  holds: 

11/  -  g\\oo  <$!=>  \\fPsi  -  gPii  ||oo  <  i. 

Let  h(u)  =|  f(u>)  —  g(u>)  |  .  Recall  that  for  a  fixed  partition,  a  single  integral  over  a 
strip,  q,  defining  the  sample  transform  takes  on  the  form 


hg  =  J2[Js  M^)  7?Mdw)  , 


where  79(u>)  is  the  characteristic  function  of  the  qth  strip.  Let  M  denote  the  area  of 
the  largest  polygon  in  our  partition,  choose  our  $i  to  be  less  than  ■£-.  That  is,  let 
the  functions  /  and  g  differ  by  no  more  than  <$i  in  the  uniform  norm  sense.  We  have 
already  proven  that  we  may  do  so  for  some  partition.  Then  we  know  for  each  element 
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of  our  sample  data  vector 

hq     =    jr([  h(u)  lq(u)du;)  , 
i=i  KJs>  J 

<  ±(—m), 

<    i.  (VI. 7) 

As  VI. 7  holds  for  each  of  the  finite  number  of  sample  integrals,  we  may  conclude  that 

\\fpS]  -qps.Woo  <i- 

Thus,  if  we  can  disregard  constraints,  VI. 4,  VI. 5,  and  VI. 6,  for  any  i  >  0,  we 
may  find  a  partition  P&x  that  ensures 

11/  "Slice       <       S1 

so  that 

||/p4l  -  gpti  ||oo    <   i- 

This  implies  that  for  any  value  of  e,  constraint  (VI. 3)  is  met. 

Temporarily  disregarding  the  constraint   \\fp  —  gp\\co  <  £5  we  have  the  less 
restrictive  optimization  problem: 

mininize:  k 

subject  to:        otj  +  k  >  g{uj)ip ^(u) ,    for  all  u  t  Q, 

— cij  +  k  >  —g(uj)ipj(u>),  for  all  well, 

clj  >  0,  for  all  j.  (VI.8) 
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That  the  value  of  the  above  optimization  problem  may  be  made  arbitrarily  small 
is  a  direct  result  of  the  fact  that  we  may  represent  any  continuous  function,  g  on 
fi  arbitrarily  well  by  a  function  of  the  desired  form,  per  Theorem  2.  Denote  the 
partition  for  which  the  value  of  k  is  less  than  e  by  Ps2.  The  term  82  is  the  largest 
largest  difference  between  two  points  in  the  same  polygon  of  P.  That  is, 

x,yesJcPs2     =>     ||s  -  y||oo    <   <^2- 

Finally,  choosing  6  to  be  the  minj^i,^},  constraint  VI. 3  is  met,  as  are  con- 
straints VI. 4  through  VI. 6  for  k  <  e.  Then  the  problem  is  feasible.  As  this  is  a  finite 
dimensional  problem,  we  employ  Theorem  13  to  ensure  that  a  solution  to  the  problem 
exists.  Therefore,  the  optimal  value  of  the  optimization  problem  is  less  than  e.        □ 

From  the  above  claim,  we  expect  that  if  our  partition  is  sufficiently  fine,  then 
we  may  reasonably  expect  to  find  an  acceptable  approximation  to  the  solution  of  our 
optimization  problem. 

C.      PROBLEM  STANDARDIZATION 

We  wish  to  understand  the  above  formulation  as  it  relates  to  our  definition  of 
the  general  linear  optimization  problem.  Before  proceeding,  it  is  vital  to  note  that 
we  are  formulating  the  problem  after  we  have  generated  a  polygonal  partition,  P,  of 
0.  Throughout  this  section,  we  assume  that  P  contains  the  n  polygons,  s\,...,sn. 
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1.        Inner  Product  Constraints:    Refining  the  Feasible 
Set 

Before  considering  the  constraints  themselves,  recall  that  the  polygonal 
partition,  P,  forms  an  rc-dimensional  basis  for  a  subset  of  the  space  of  functions  from 
which  we  select  our  optimal  function.  We  may,  consequently,  think  of  any  density 
function  /  as  a  vector  y  c  7£n,  where  the  jth  component  of  y  is  the  scalar  value  of  the 
density  on  polygon,  Sj.  For  reasons  that  become  clear  shortly,  we  augment  the  vector 
of  decision  variables  to  be  y  =  [qi,  q2,  . . .  ,an,  k]T  e  7ln+l . 

We  divide  our  constraints  into  three  distinct  classes: 

1)  Strip  based  constraints, 

2)  Polygon  based  constraints,  and 

3)  fl  based  constraints. 

First  consider  the  strip  based  constraints.  We  require  that  the  sample 
transform  of  the  optimal  objective  density  be  within  a  specified  tolerance  of  the 
observed  sample  transform.  The  constraints  were  identified  in  the  previous  section 
by  the  equation 

ll/p-blloo     <    £. 

Eliminating  the  norm  above  results  in  the  two  constraints 

-fp    >    -b-e,  (VI.9) 

and 

fp     >    b-£.  (VI. 10) 
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Let  m  denote  the  number  of  strips  used  to  generate  the  partition.  P. 
Let  Q  —  {<7i,<72, . . .  ,qm}  be  the  set  of  such  strips.  Then  for  each  q  c  Q,  we  require 
that  the  jth  component  of  our  constraint  vector  be  determined  by  the  following  rule: 

afai)     =    area(sj)    max{7q(w)  V>j(w)}  . 

Of  course,  this  convention  is  the  same  as  that  of  Chapter  III.  The  jth  component  of 
the  vector,  a(q;)  is  the  area  of  the  jth  polygon  if  the  polygon  falls  within  strip  qt,  and 
zero  otherwise.  As  before,  we  wish  to  consider  the  two  constraints  associated  with 
each  strip  separately,  and  define  Q~  and  Q+  to  index  constraints  (VI. 9)  and  (VI.  10) 
respectively. 

The  right  hand  side  of  each  constraint  is  also  determined  as  in  Chapter 
III.  That  is,  b(q— )  and  b(q+)  as  in  equations  (VI. 9)  and  (VI. 10),  where  6,  is  the  data 
from  strip  g,  of  our  sample  transform.  We  append  a  zero  to  each  strip  based  constraint 
vector,  as  each  is  independent  of  the  value  k. 

The  polygon  based  constraints  are  found  entirely  in  the  requirement 
that  our  density  function  be  non-negative.  In  the  initial  formulation,  the  requirement 
was  written 

otj     >     0,       for  all  j. 

That  is,  we  require  that  the  density  assigned  to  each  polygon  in  the  optimal  vector 
be  non-negative.  Let  P  =  {si,  52, ... ,  sn}  be  the  fixed  partition.  Then  the  constraint 
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vector  associated  with  each  polygon  is  formed  by  the  simple  rule: 

aj(si)     =     max{if>j(u))  il>i(u)}, 

J   1,     if  *  =  j 

0,      otherwise, 

for  ij  =  l,2,...,n. 

We  append  a  zero  to  each  a(si)  as  in  the  case  of  the  strip  based  con- 
straints, as  we  are  selecting  a  vector  from  7Zn+i .  Clearly,  6(5,)  =  0,  for  all  i.  Then  the 
vector  form  of  each  polygon  based  constraint  is: 

(a(si),y)   >  0,  for  i  =  l,2,...,n. 

As  our  third  class  of  constraints  corresponds  to  the  set  $7,  we  may 
correctly  infer  that  the  final  index  sets  are  infinite.  These  index  sets  provide  the 
constraints  that  facilitate  a  comparison  of  solution  quality.  There  are,  in  fact,  two 
such  index  sets,  $7+  and  fi~,  as  we  again  eliminate  the  explicit  use  of  the  infinity 
norm  from  the  formulation.  In  the  initial  problem  statement,  these  constraints  were 
written 


Qj  +  k     >    g(u)  t/>j(u>), 
— a j  +  k     >        —g{u)  ijjj(uj),       for  all  u  t  Cl. 
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We  focus  only  on  the  former,  as  Q~  is  formulated  in  a  nearly  identical 
manner,  and  the  process  has  been  executed  in  the  strip  based  constraints.  We  desire 
constraints  of  the  form  (a(cj+),y)    >  b(u;+),  for  all  u>+  e  Q+ .    Let 


a,(u,+  )     =  ^{u+),      forj  =  l,2,...,n,  u+  e  ft+, 
an+i(u;/")  =  1,  for  all  u>+  e  ft+ . 

The  associated  b(u>+)  is  defined  to  be  g(u+).  The  constraints  associated  with  the  set 
fi~  are  formed  in  exactly  the  same  manner,  with  sign  changes  as  appropriate. 

Concluding,  we  define  our  index  set  T  =  Q+  U  Q~   U  P  U  n+  U  Q~. 

D.      THE  IMAGE  RECONSTRUCTION  DUAL 

The  image  reconstruction  optimization  problem,  as  we  have  formed  it,  is  a 
specific  example  of  the  uniform  approximation  problem.  Consequently,  we  find  some 
strong  similarities  in  its  dual  problem  to  the  dual  of  the  approximation  of  the  ex- 
ponential function.  Let  us  derive  the  dual,  Z),  of  our  image  reconstruction  problem. 
Note  that  this  section  is  included  in  the  interest  of  completeness.  The  material  herein 
is  complicated  and  is  not  especially  enlightening.  The  reader  may  wish  to  skip  this 
section. 

We  seek  a  subset,  {^,  t2, . . .  ,tg}  C  T,  and  the  non-negative  vector  x  = 
[xi ,  X2, . .  • ,  xq]     that  maximize  the  equation 

(b,x) 
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and  satisfies 

(x,a(t))     =     c. 

We  address  the  selection  of  the  subset  first.  Consider  the  strip-based  con- 
straints, associated  with  Q  C  T.  Recall  that  Q  =  Q+  U  Q~ .  We  seek  some  subset  of 
each  of  these  sets.  We  denote  these  subsets  Q+,  and  Q~ .  With  each  element  of  each 
subset,  we  associate  a  non-negative  real  number,  x(q+),  for  j  =  1, . . . ,  ng-+ ,  and  x(qj), 
for  j  =  l,...,n,-. 

Considering  the  index  set,  P,  with  which  we  associated  the  polygon  based 
constraints,  we  seek  some  non-negative  value  x(sj)  to  associate  with  each  constraint 
of  a  subset  P  C  P.  Following  the  above  convention,  we  let  j '  =  1, ... ,  np. 

Let  us  move  to  the  infinite  index  sets,  f2+  and  fi~ '.  As  we  noted  above,  there 
are  two  classes  of  constraint  vectors  associated  with  our  index  sets,  Q,+  and  fi~.  In 


particular, 


1 


and 


-ipi{u    ) 
— 02(u>") 

1 


For  each  of  our  index  sets,  Q+  and  Cl   ,  we  seek  some  discretization 

fi+     =     {u>+,u;+,...,u;+ft+},  and 

Cl~     =     {wf  ,a;J, . . .  ,u;~  _},  as  well  as  non-negative  scalars, 
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x{wf),x(u}),...,x(u+^)t  and 

Then  the  dual  D  is  to  find  the  above  discretization  and  non-negative  x  values 
that  maximize  the  expression: 

E  Aqtmt)  -  E  *(*  w«r)  +  E  *(*)*(«)  +  E  *(»t)K<*t)  -  E  *(*rW"r), 

t=i  i=i  1=1  i=i  i=i 

while  satisfying  the  constraints: 

nQ+  nQ-  nf>  nn+  nn- 

E  x^t)ar{qt  )-E  *(9f  )Mgr)+E*(s«X(s.-)+E  ^(^K^  )~e  *(wf")ar(wr)  =  o, 

i=i  1=1  i=i  t=i  t=i 

for  r  =  1,2, . . .  ,  n,  and 


"<5+  nQ- 

E  ^(^K+iUf )  -  E  *(%"  K+ifaf ) 

1=1  t=i 


+  EX(5')an+l(5t)  +   E  X(Ut)an+l(u?) 


t=l  i=l 

i=l 

z(9,+)  >  0,      for     i  =  l,2,...,ri0+, 
*(?,")  ^  °>     for     >  =  1,2,  ...,ng_, 
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x(si)  >  0,      for      i  =  1,2,. . .  ,rap, 

^(^,+  )  >  05      for      i=l,2 n^+,  and. 

s(a>f)  >  0,      for      i  =  1,2, . . .  ,n^_. 

While  the  above  formulation  of  the  dual  is  intimidating,  we  may  simplify  im- 
mediately by  recognizing  some  features  of  the  constraints  of  our  primal  problem.  We 
know  that  the  scalar,  6(s,)  =  0,  for  all  st.  Then  the  middle  term  in  the  objective 
function  disappears  completely. 

Let  us  move  to  the  first  constraint.  The  middle  sum  also  collapses  to  the  single 
term  x(sr),  as  we  have  defined  ar(sj)  to  be  the  Kronecker  6(r,j).  The  first  three 
terms  of  the  second  constraint  disappear  altogether,  as  we  have  specified,  an+\(sj)  = 
an+i(sj)  =  0,  for  all  j.  The  non-negativity  constraints  remain  the  same. 

E.      A  SAMPLE  SOLUTION 

We  now  use  the  formulation  of  the  image  reconstruction  problem  as  a  linear 
program  to  solve  a  simple  problem.  We  first  discuss  the  geometry  of  the  partition 
that  we  are  using,  and  then  identify  some  additional  simplifying  assumptions  that 
make  the  problem  more  tractable.  We  introduce  the  expected  density  of  our  sample 
problem,  and  conclude  with  the  Simplex  solution  of  the  problem. 

1.       The  Partition 

The  partition  that  we  use  in  this  example  is  illustrated  in  Figure  16, 
where  the  color  of  a  polygon  is  a  function  of  its  area.    Larger  values  correspond  to 


Figure  16.    The  Partition  of  the  Sample  Problem 


lighter  colors.  We  have  chosen  the  four  angles  0,  ^,  |,  and  ^L,  with  five  strips  for  each 
angle.  The  resulting  partition  consists  of  89  polygons.1  It  should  be  noted  that  each 
strip  has  width  of  ~.  Consequently,  as  views  at  angles  of  ^  and  —■  require  more  than 
5  strips  to  cover  the  unit  square  completely,  only  the  portion  of  the  square  that  falls 
in  the  five  center  strips  is  considered.  The  rest  is  ommitted. 

2.       A  Simplifying  Assumption 

Rather  than  attempt  to  solve  the  infinite  dimensional  problem  as  de- 
rived in  the  initial  portion  of  this  chapter,  we  project  the  target  density  onto  the 
n-dimensional  polygonal  basis  of  our  partition.    That  is,  we  insist  that  the  target 


JThe  manner  in  which  the  polygons  were  identified  and  the  areas  of  each  polygon  computed  is 
not  a  matter  of  particular  concern  here.  It  is  sufficient  to  state  that  the  symmetry  of  the  partition 
was  deeply  exploited  in  a  manner  which  simplifies  the  problem  of  polygon  identification  and  area 
computation  over  89  polygons  to  one  of  many  fewer  than  89. 
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function  be  constant  on  each  of  the  polygons  of  the  partition.  This  simplification 
reduces  the  infinite  index  sets  0+  and  fi~  to  finite  sets,  as  we  need  only  consider  a 
representative  ujj  c  Sj  for  each  polygon  s^  when  determining  the  norm  of  the  differ- 
ence between  our  target  function  and  optimal  function.  Without  this  assumption,  the 
problem  is  very  similar  to  the  infinite  problem  of  approximating  the  exponential  with 
polynomials,  which  was  discussed  in  more  detail  when  the  Simplex  algorithm  was 
introduced.  It  is  possible  that  this  problem  is  solveable  without  this  simplification, 
but  no  attempt  is  made  to  solve  it  in  this  thesis. 

We  choose,  in  projecting  the  target  function  onto  our  finite  dimensional 
space,  the  density  of  the  function  over  each  polygon  divided  by  the  area  of  the  polygon. 
That  is,  after  the  target  function  g  is  projected  into  the  finite  space,  it  takes  on  the 
form: 

n 

g  =  ^Pjifrj, 
where 

That  is,  )3  =  the  mass  of  g  over  polygon  Sj,  divided  by  the  area  of  polygon  Sj. 

3.       The  Target  Function 

We  now  identify  the  target  function  of  the  sample  problem.  We  use  the 
simple  function, 

g(x,y)=  ^ar  —  -J    +  (y  -  -J    , 
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Figure  17.    The  Projected  Target  Function 


for  x  t  [0,1]  and  y  e  [0,1].  The  projection  of  the  density  function  is  illustrated  in 
Figure  17.  The  particular  data  for  the  constant  densities  assigned  to  each  of  the 
strips  represent  the  values  which  we  hope,  or  expect,  to  find  in  the  solution  of  our 
problem,  before  considering  the  data.  That  is,  the  values  are  assumed  to  represent 
the  most  likely  solution  to  our  problem. 

4.       The  Test  Density 

The  density  function  that  we  use  to  generate  the  test  data  is  given  by 


the  expression 


h{x,y)    = 


{*-\)2  +  {y-\)\   f<*  (*-^)2<A 


1,  otherwise. 

The  density  function  is  displayed  in  Figure  18.     The  values  of  sample  transform 
defining  integrals  become  the  right  hand  side  of  our  equality  constraints  when  we 
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Figure  18.   The  Test  Density  Function 


formulate  the  problem. 

The  manner  in  which  we  have  defined  the  projection  of  a  density  onto 
the  finite  dimesional  space  assures  us  that  both  the  test  density  h  and  its  projection 
have  the  same  sampled  transform.  Thus,  barring  catastrophic  rounding  error,  the 
formulation  is  always  feasible.  That  is,  there  must  be  some  density  function  that 
produces  the  sampled  transform,  even  after  we  have  projected  the  test  density  onto 
the  partition.  If  the  sampled  transform  is  uniquely  determined,  we  reconstruct  the 
projection  perfectly,  though  the  value  of  the  variable  k  may  be  quite  large. 

As  a  basis  of  comparison,  we  note  that  the  maximum  difference  between 
the  projection  of  the  target  density  and  the  test  density  is  given  by  d  =  .1949.  We 
may  certainly  expect  then,  that  the  optimal  density  varies  by  no  more  than  the  above 
value  of  d. 
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Figure  19.    The  Simplex  Solution  of  the  Sample  Problem 


The  optimal  density  as  determined  by  the  Simplex  algoritm  is  displayed 
in  Figure  19.  The  value  achieved  for  the  maximum  absolute  deviation  between  the 
target  density  and  the  optimal  density  is  d'  =  .1577.  We  consider  the  difference  over 
each  polygon  in  Figure  20. 


F.       SUITABILITY   OF   SIMPLEX   IN   IMAGE   RECON- 
STRUCTION 

We  briefly  consider  the  merit  of  using  the  Simplex  algorithm  to  solve  the  image 
reconstruction  problem.  That  is,  we  wish  to  consider  how  well  the  tool  we  have  chosen 
fits  our  particular  job. 

The  results  of  this  particular  example  show  the  tendency  of  this  formulation 
to  spread  error  over  the  entire  region.  This  consequence,  it  is  believed,  results  from 
the  use  of  the  infinity  norm.    We  may  also  question  the  choice  of  target  functions, 
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Figure  20.  Difference  over  Each  Polygon 


and  may  look  to  other  methods  of  qualification.  However,  the  optimization  problem 
achieves  what  it  is  designed  to  achieve.  That  is,  we  have  found  the  density  that 
satisfies  the  minimum  deviation  in  the  uniform  norm  sense. 

We  are  also  forced  to  consider  the  substantial  data  that  are  required  to  solve  the 
problem.  The  problem  of  polygon  identification  is  a  difficult  one  by  itself,  especially 
in  view  of  the  fact  that  a  typical  partition  for  the  CAT  scan  problem  is  generated  by 
200-300  angles  with  up  to  500  strips  per  angle.  With  this  geometry,  we  know  that  the 
number  of  polygons  exceeds  4,000,000.  Further,  we  require  the  area  of  each  polygon 
be  known  to  solve  the  problem  as  we  have  formulated  it.  Finally,  as  each  polygon 
gives  rise  to  a  variable  in  our  primal  problem,  we  are  solving  a  problem  in  a  subspace 
of  7Zn  where  n  is  quite  large.  On  the  positive  side,  we  know  that  we  must  solve  the 
polygon  identification  and  polygonal  area  problems  only  once.   Further,  the  matrix, 
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A  which  results  from  the  formulation  above  is  extremely  sparse,  which  may  lead  to  a 
more  rapid  solution  of  the  Simplex  problem,  or  invite  other  methods  of  solving  Linear 
Programs. 

In  conclusion,  the  author  contends  that  the  Simplex  algorithm  fits  well  con- 
ceptually, but  may  not  suited  for  the  vastness  of  the  problem  as  it  is  formulated  here. 
Projecting  the  density  functions  onto  the  polygonal  partion  is  conceptually  identical 
to  selecting  finite  subsets  of  an  infinite  index  set.  The  theorems  presented  in  regard  to 
the  image  reconstruction  problem  indicate  that  we  may  solve  the  infinite  dimensional 
problem  through  a  sequence  of  finite  dimensional  problems,  when  certain  conditions 
are  met. 

Some  alternatives  that  might  warrant  future  consideration  are  norms  other 
than  the  infinity  norm,  or  using  the  Simplex  model  to  refine  existing  solutions  to  the 
reimaging  problem,  where  the  number  of  variables  is  less  restrictive. 
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VII.         CONCLUSION 

In  conclusion,  we  have  introduced  the  Simplex  algorithm  in  a  context  quite 
apart  from  its  usual  applications.  The  principal  vehicles  for  the  exploration  of  the 
algorithm  were  three  unique  applications,  each  illuminating  distinct  features  of  the 
theory  underlying  the  implementation  of  the  Simplex  algorithm. 

In  particular,  we  began  with  a  problem  of  finding  orthogonal  monic  polynomi- 
als over  a  closed  interval.  This  example  led  to  a  very  basic  Simplex  formulation,  and 
was  solved  as  a  finite  dimensional  problem.  The  requirement  that  the  polynomials 
be  monic  facilitated  the  relatively  simple  formulation.  Follow  on  problems  to  this 
example  might  be  the  adaptation  of  the  algorithm  to  generate  an  non-polynomial 
orthogonal  basis  for  an  infinite  dimensional  function  space,  or  perhaps  to  fit  the  al- 
gorithm to  solving  the  non-linear  orthonormal  basis  generation  problem. 

Second,  we  formulated  the  problem  of  approximating  a  function  over  a  closed 
interval  in  the  uniform  norm  sense.  Unlike  the  first  example,  the  problem  was  infinite 
dimensional,  in  that  the  formuation  required  a  constraint  for  every  number  in  the 
uncountably  infinite  set.  This  problem  proved  particularly  helpful  in  illustrating  the 
principle  of  weak  duality,  and  ultimately  illustrated  the  Simplex  algorithm  itself.  The 
special  qualities  of  polynomial  approximation  were  ommitted,  though  the  reader  is 
referred  to  [Ref.  9]  for  a  more  complete  discussion  thereof.  Again,  potential  areas  for 
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future  research  might  include  approximation  with  functions  other  than  polynomials. 

Each  of  the  above  examples  were  used  extensively  to  illustrate  the  highlights  of 
convexity  and  duality,  upon  which  the  Simplex  algorithm  is  based.  The  treatment  was 
relatively  general,  though  many  of  the  theorems  required  that  the  linear  optimization 
problem  be  finite.  Work  is  underway  to  identify  classes  of  infinite  dimensional  prob- 
lems which  may  be  solved  by  a  sequence  of  finite  dimensional  problems.  The  reader 
is  referred  to  [Ref.  11]  for  more  complete  discussion  of  this  active  area  of  research. 
Highlights  include  infinite  horizon  planning,  fuzzy  set  semi-infinite  programming,  and 
linear  programming  in  control  theory. 

Another  area  of  focus  in  this  paper  was  on  the  Image  Reconstruction  problem. 
Again,  this  is  an  area  of  active  research.  After  presenting  the  requisite  background, 
we  formulated  this  problem  as  an  infinite  dimensional  linear  optimization  problem, 
and  as  a  special  case  of  the  uniform  approximation  problem.  Results  were  presented 
that  indicated  that  use  of  Simplex  to  solve  a  sequence  of  linear  programs  is  conceptu- 
ally sound,  though  not  necessarily  practical  in  view  of  the  size  of  the  problem.  This 
application  of  the  algorithm,  however,  is  open  to  more  extensive  research  in  a  number 
of  areas.  A  different  choice  of  norms  by  which  the  quality  of  density  functions  is 
measured  may  eliminate  a  number  of  constraints.  A  technique  for  formulating  opti- 
mization problems  with  the  2-norm  is  found  in  [Ref.  12],  and  may  prove  useful  in 
this  application.  The  Simplex  algorithm  may  also  provide  an  inexpensive  method  to 
solve  coarser  problems,  from  which  one  may  determine  the  necessity  of  constructing 
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more  detailed  models.  Alternately,  there  may  be  some  utility  in  using  the  algorithm 
to  solve  the  reconstruction  problem  in  only  small  portions  of  the  set  over  which  a 
density  is  defined.  If  there  is  utility  in  such  an  application,  the  logical  consequence  is 
research  of  parallel  Simplex  implementation. 

The  potential  utility  of  the  Simplex  algorithm  to  unconventional  applications 
seems  clear.  Even  when  actual  implementation  of  the  algorithm  is  not  practical,  the 
tools  of  convexity  and  duality  apply  to  broader  areas  of  optimization. 
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